### Arguments ###
## shift.args: a list object to be specified to simulate and test clade-associated diversification. 
## If NULL or empty list() the test is not performed. The list must include:
  # BMrates = set as TRUE to perform search.shiftDR by computing phenotypic evolutionary rates by using geomorph's compare.evol.rates
  # RRrates = set as TRUE to perform search.shiftDR by computing phenotypic evolutionary rates by using RRphylo
  # auto = set as TRUE to automatically scan the phylogeny looking for clades showing significant rate shifts (in both DR and phenotypic evolutionary rates)
  # node = set as TRUE to test the clades imposed to shift
  # overfit = to test the robustness of search.shiftDR results to sampling by applying the funcion overfitDR. 
  #           To perform overfit for search.shiftDR as applied either under BMrates or under RRrates or both, 
  #           the user needs to specify the arguments clus and nsim as follows: 
  #           overfit=list(BMrates=list(nsim=10,clus=0.5),RRrates=list(nsim=10,clus=0.5))
##
## trend.args: a list object to be specified to simulate and test character-associated diversification. 
## If NULL the test is not performed. To perform search.shiftDR only (i.e. not testing the effect of sampling),
## the list must be empty list(). To test the robustness of search.shiftDR results to sampling by applying 
## the funcion overfitDR, the user must specify the arguments clus and nsim as follows: 
## trend.args=list(overfit=list(nsim=10,clus=0.5)). ESsim = TRUE to compare search.shiftDR to ES-sim
##
## reps = the number of iterations

### Results ###
## The function returns summary results tables.
## If clade-associated diversification is tested (i.e. shift.args is not NULL) under the automatic mode,
## the table includes (for either BMrates=TRUE or RRrates=TRUE)
  # empty =  the percantege of simulations producing no significant results.
  # double = the percentage of simulations where another clade in addition to the clade simulated to shift was found. 
  #          These instances are also counted as significant results.
  # error = the percentage of simulations where the function found the wrong clade to shift (not the one we simulated to shift).
  # DRpBRp-DRpRRp = the percentage of simulations where the clade simulated to shift was found to have positive shift in both DR and BR/RR rates.
  # DRpBRm-DRpRRm = the percentage of simulations where the clade simulated to shift was found to have positive shift in DR andin BR/RR rates.
  # DRmBRp-DRmRRp = the percentage of simulations where the clade simulated to shift was found to have negative shift in DR and positive shift BR/RR rates.
  # DRmBRm-DRmRRm = the percentage of simulations where the clade simulated to shift was found to have negative shift in both DR and BR/RR rates.
  # type1 = the type 1 error rate
##
## If clade-associated diversification is tested (i.e. shift.args is not NULL) by selecting a clade to be analyzed,
## the output includes (for either BMrates=TRUE or RRrates=TRUE) two result tables. In the $focal table, results for the focal clade are summarized.
## For simulations where focal clades produced no significant results, the node along the path aren inspected ($path). The tables include: 
  # error = the percentage of simulations  producing no significant results.
  # DRpBRp-DRpRRp = the percentage of simulations where the clade simulated to shift was found to have positive shift in both DR and BR/RR rates.
  # DRpBRm-DRpRRm = the percentage of simulations where the clade simulated to shift was found to have positive shift in DR andin BR/RR rates.
  # DRmBRp-DRmRRp = the percentage of simulations where the clade simulated to shift was found to have negative shift in DR and positive shift BR/RR rates.
  # DRmBRm-DRmRRm = the percentage of simulations where the clade simulated to shift was found to have negative shift in both DR and BR/RR rates.
  # type1 = the type 1 error rate
##
## If character-associated diversification is tested (i.e. trend.args is not NULL), the table includes:
  # power =  the percentage of simulations where a positive (negative) association between DR and phenotype was imposed 
  #          which produce significant and positive (negative) slopes.
  # misclass =  the percentage of simulations where a positive (negative) association between DR and phenotype was imposed 
  #          which produce significant and negative (positive) slopes.
  # type1 = the type 1 error rate

## Lists of results will be stored into the global environment depending on the function settings.
## If clade-associated diversification is tested (i.e. shift.args is not NULL) 
## the following elements might be stored into the global environment, depending on settings 
## (i.e. wheter BMrates, RRrates, auto, and node are TRUE and overfit is not NULL):
  # tree.list = list of phylogenetic trees generated to simulate a positive shift in DR
  # tree2.list = list of phylogenetic trees generated to simulate a negative shift in DR
  # y.list = list of phenotypic values generated to simulate a positive shift in DR (related to tree.list)
  # y2.list = list of phenotypic values generated to simulate a negative shift in DR (related to tree2.list)
  # yB.list = list of phenotypic values generated to assess the Type 1 error rate on the tree generated to simulate a positive shift in DR (related to tree.list)
  # yBm.list = list of phenotypic values generated to assess the Type 1 error rate on the tree generated to simulate a negative shift in DR (related to tree2.list)
  # NN = node subtending to the clade set to show a positive shift in DR rate (related to tree.list)
  # NN2 = node subtending to the clade set to show a negative shift in DR rate (related to tree2.list)
  # sigmas = the phenotypic evolutionary rate value imposed at the shifting clade
  # DRratio = the ratio between the mean DR of species within the clade set to show a positive DR shift before 
  #           modification and the mean DR of the rest of the tree, divided by the same figure after clade modification (related to tree.list)
  # DRratio.m = the ratio between the mean DR of species within the clade set to show a negative DR shift before 
  #             modification and the mean DR of the rest of the tree, divided by the same figure after clade modification (related to tree2.list)
  # PH = the proportion of tree tips to be added or removed to test the effect of sampling and anagenesis within overfitDR
  # res.list-resRR.list-res.node.list-resRR.node.list = summary results list for BMrates automatic, RRrates automatic, BMrates node, and RRrates node, respectively.
  # SSdr.list-SSdrB.list-SSdrm.list-SSdrBm.list = results for search.shiftDR under BMrates automatic for test of positive DR simulations, 
  #                                               test of type 1 error on positive DR simulations, test of negative DR simulations, 
  #                                               test of type 1 error on negative DR simulations, respectively.
  # SSdrRR.list-SSdrBRR.list-SSdrmRR.list-SSdrBmRR.list = results for search.shiftDR under RRrates automatic for test of positive DR simulations, 
  #                                                       test of type 1 error on positive DR simulations, test of negative DR simulations, 
  #                                                       test of type 1 error on negative DR simulations, respectively. 
  # SSdr.node.list-SSdrB.node.list-SSdrm.node.list-SSdrBm.node.list = results for search.shiftDR under BMrates node for test of positive DR simulations, 
  #                                                                   test of type 1 error on positive DR simulations, test of negative DR simulations, 
  #                                                                   test of type 1 error on negative DR simulations, respectively.
  # SSdrRR.node.list-SSdrBRR.node.list-SSdrmRR.node.list-SSdrBmRR.node.list = results for search.shiftDR under RRrates node for test of positive DR simulations, 
  #                                                                           test of type 1 error on positive DR simulations, test of negative DR simulations, 
  #                                                                           test of type 1 error on negative DR simulations, respectively.
  # OFdr.list-OFdrB.list-OFdrm.list-OFdrBm.list = results for overfitDR under BMrates automatic for test of positive DR simulations, 
  #                                               test of type 1 error on positive DR simulations, test of negative DR simulations, 
  #                                               test of type 1 error on negative DR simulations, respectively.
  # OFdrRR.list-OFdrBRR.list-OFdrmRR.list-OFdrBmRR.list = results for overfitDR under RRrates automatic for test of positive DR simulations, 
  #                                                       test of type 1 error on positive DR simulations, test of negative DR simulations, 
  #                                                       test of type 1 error on negative DR simulations, respectively.
  # OFdr.node.list-OFdrB.node.list-OFdrm.node.list-OFdrBm.node.list = results for overfitDR under BMrates node for test of positive DR simulations, 
  #                                                                   test of type 1 error on positive DR simulations, test of negative DR simulations, 
  #                                                                   test of type 1 error on negative DR simulations, respectively.
  # OFdrRR.node.list-OFdrBRR.node.list-OFdrmRR.node.list-OFdrBmRR.node.list = results for overfitDR under RRrates node for test of positive DR simulations, 
  #                                                                           test of type 1 error on positive DR simulations, test of negative DR simulations, 
  #                                                                           test of type 1 error on negative DR simulations, respectively.
##
## If character-associated diversification is tested (i.e. trend.args is not NULL) 
## the following elements might be stored into the global environment, depending on settings:
  # tree.trend.list = list of phylogenetic trees generated to test for character associated diversification.
  # y.trend.list = list of phenotypic values generated to test for character associated diversification.
  # yB.trend.list = list of phenotypic values generated to assess the type 1 error rate in character associated diversification.
  # res.trend.list = summary results list for search.shiftDR.
  # PH.trend = the proportion of tree tips to be added or removed to test the effect of sampling and anagenesis within overfitDR.
  # Div.trend = proportion of tree species sampled to be added a sister species.
  # Esp.trend = exponent to condition the sampling probability.
  # tree.ratio.trend = ratio of the number of tree tips before and after modifications.
  # means.shift.trend =  percentile of the mean of the final phenotypic vector relative to the original one.
  # resOF.trend.list = summary results list for overfitDR.
  # SSdr.trend.list-SSdrB.trend.list = results of search.shiftDR for character-associated
  #                                     diversification test and type 1 error test, respectively.
  # OFdr.trend.list-OFdrB.trend.list = results of overfitDR for character-associated 
  #                                    diversification test and type 1 error test, respectively.


required<-c("RRphylo","ape","geiger","phytools","picante","stats4","R.utils","Rphylopars",
            "scales","geomorph","plotrix","RColorBrewer","doParallel","parallel","graphics",
            "mvtnorm","phangorn")

if(any(!required%in%installed.packages()[,1])){
  install.packages(required[which(!required%in%installed.packages()[,1])])
}

diversification.test<-function(shift.args=list(BMrates=TRUE,RRrates=TRUE,auto=TRUE,node=FALSE,overfit=list(BMrates=NULL,RRrates=NULL)),
                     trend.args=list(overfit=NULL,ESsim=FALSE),reps=10){
  
  require(RRphylo)
  require(ape)
  require(geiger)
  require(phytools)
  require(picante)
  require(stats4)
  require(R.utils)
  require(Rphylopars)
  require(mvtnorm)
  require(phangorn)
  
  if(!is.null(shift.args)){
    tree.list<-tree2.list<-y.list<-y2.list<-yB.list<-yBm.list<-list()
    NN<-NN2<-sigmas<-DRratio<-DRratio.m<-PH<-array()
    res.list<-resRR.list<-res.node.list<-resRR.node.list<-list()
    SSdr.list<-SSdrB.list<-SSdrm.list<-SSdrBm.list<-list()
    SSdrRR.list<-SSdrBRR.list<-SSdrmRR.list<-SSdrBmRR.list<-list()
    SSdr.node.list<-SSdrB.node.list<-SSdrm.node.list<-SSdrBm.node.list<-list()
    SSdrRR.node.list<-SSdrBRR.node.list<-SSdrmRR.node.list<-SSdrBmRR.node.list<-list()
    OFdr.list<-OFdrB.list<-OFdrm.list<-OFdrBm.list<-list()
    OFdrRR.list<-OFdrBRR.list<-OFdrmRR.list<-OFdrBmRR.list<-list()
    OFdr.node.list<-OFdrB.node.list<-OFdrm.node.list<-OFdrBm.node.list<-list()
    OFdrRR.node.list<-OFdrBRR.node.list<-OFdrmRR.node.list<-OFdrBmRR.node.list<-list()
  }
  
  if(!is.null(trend.args)){
    ecdf_fun<-function(x,perc) ecdf(x)(perc)
    
    tree.trend.list<-y.trend.list<-yB.trend.list<-
      res.trend.list<-resOF.trend.list<-list()
    PH.trend<-Div.trend<-Esp.trend<-tree.ratio.trend<-means.shift.trend<-array()
    SSdr.trend.list<-SSdrB.trend.list<-list()
    OFdr.trend.list<-OFdrB.trend.list<-list()
    ESsims.list<-ESsimsB.list<-res.ESsim.list<-list()
  }
  
  i=1
  while(i<=reps){
    print(i)
    ##### Clade-associated diversification #####
    if(!is.null(shift.args)){
      
      shift.args$BMrates->BMrates
      shift.args$RRrates->RRrates
      shift.args$auto->auto
      shift.args$node->node
      shift.args$overfit->Soverfit
      
      if(is.null(BMrates)) BMrates<-TRUE
      if(is.null(RRrates)) RRrates<-TRUE
      if(is.null(auto)) auto<-TRUE
      if(is.null(node)) node<-FALSE
      
      repeat({
        try({  
          ### creating a phylogenetic tree with > 180 species
          repeat({
            sim.bdtree(b = .5, d = 0.2)->tree
            if(Ntip(tree)>180) break
          })
          
          ### computing diversifiction rates for the tree
          as.matrix(1/evol.distinct(tree)[,2])->DRt
          rownames(DRt)<-evol.distinct(tree)[,1]
          
          #### Procedure to generate a positive shift in DR ####
          ### selecting a clade to be replaced
          sample(seq(1/10,1/5,length.out = 10),1)->Cdiv
          round(Ntip(tree)*Cdiv)->Csize
          subtrees(tree)->subtr
          sapply(subtr,function(x) getMRCA(tree,x$tip.label))->names(subtr)
          if(any(sapply(subtr,Ntip)==Csize))
            subtr[which(sapply(subtr,Ntip)==Csize)]->Ctar else
              subtr[which(sapply(subtr,Ntip)>Csize*0.9&sapply(subtr,Ntip)<Csize*1.1)]->Ctar
          if(length(Ctar)<1) subtr[which.min(abs(Csize-sapply(subtr,Ntip)))]->Ctar
          
          if(length(Ctar)>1) Ctar[sample(length(Ctar),1)]->Ctar
          Ctar[[1]]->Ctar
          getMRCA(tree,Ctar$tip.label)->nodo
          tree$edge.length[which(tree$edge[,2]==nodo)]->rel
          max(nodeHeights(Ctar))->Htar
          
          ### generating a new tree being twice the size of the selected clade at least
          withTimeout(repeat({
            sim.bdtree(b =.8, d = 0.2)->Cnew
            if(Ntip(Cnew)>=2*Ntip(Ctar) & Ntip(Cnew)<.5*(Ntip(tree)-Ntip(Ctar)+Ntip(Cnew))) break
          }),timeout=60)
          
          ### rescaling the newly generated tree on the height of the selected clade
          ### and binding it to the phylogeny in place of such clade
          geiger::rescale(Cnew,"depth",Htar)->Cnew
          
          gsub("s","t",Cnew$tip.label)->Cnew$tip.label
          Cnew$root.edge<-0
          bind.tree(tree,Cnew,nodo)->temp
          drop.tip(temp,Ctar$tip.label)->tree1->tree.list[[i]]
          getMRCA(tree1,Cnew$tip.label)->N->NN[i]
          
          ### simulating a phenotype according to the Brownian motion for the whole tree
          fastBM(tree1)->y1
          fastAnc(tree1,y1)[as.character(N)]->a
          
          ### simulating a phenotype for the binded clade
          ### by sampling the sig2 in between 0.1-0.3 and 3-10, so that it is different
          ### from the rate of the rest of the tree
          sigg<-sample(c(seq(.1,.3,length.out = 100),seq(3,10,length.out = 100)),1)
          sigg->sigmas[i]
          fastBM(Cnew,a=a,sig2=sigg)->ynew
          y1[tips(tree1,N)]<-ynew
          y1->y.list[[i]]
          
          ### computing diversifiction rates for the whole tree and comparing them to
          ### the DR of the original tree
          as.matrix(1/evol.distinct(tree1)[,2])->DRs
          rownames(DRs)<-evol.distinct(tree1)[,1]
          
          mean(DRt[match(Ctar$tip.label,rownames(DRt))])/mean(DRt[-match(Ctar$tip.label,rownames(DRt))])->rtp
          mean(DRs[match(Cnew$tip.label,rownames(DRs))])/mean(DRs[-match(Cnew$tip.label,rownames(DRs))])->rs
          rs/rtp->DRratio[i]
          
          #### Repeating the procedure on the original tree to generate a negative shift in DR ####
          ### selecting a clade to be replaced
          names(sample(which(sapply(subtr,Ntip)>40 & sapply(subtr,Ntip)<.55*Ntip(tree)),1))->Ctar2
          extract.clade(tree,as.numeric(Ctar2))->Ctar.m
          max(nodeHeights(Ctar.m))->Htar.m
          
          ### generating a new tree being 25% the size of the selected clade at most
          withTimeout(repeat({
            sim.bdtree(b =.5, d = 0.1,stop="taxa",n=.25*Ntip(Ctar.m))->Cnew.m
            if(Ntip(Cnew.m)>=.25*Ntip(Ctar.m)) drop.tip(Cnew.m,sample(Cnew.m$tip.label,(Ntip(Cnew.m)-floor(.25*Ntip(Ctar.m)))))->Cnew.m
            if(Ntip(Cnew.m)>=10) break
          }),timeout=60)
          as.numeric(Ctar2)->nodo2
          
          ### rescaling the newly generated tree on the height of the selected clade
          ### and binding it to the phylogeny in place of such clade
          geiger::rescale(Cnew.m,"depth",Htar.m)->Cnew.m
          gsub("s","t",Cnew.m$tip.label)->Cnew.m$tip.label
          bind.tree(tree,Cnew.m,nodo2)->temp
          drop.tip(temp,Ctar.m$tip.label)->tree2->tree2.list[[i]]
          getMRCA(tree2,Cnew.m$tip.label)->N2->NN2[i]
          
          ### simulating a phenotype according to the Brownian motion for the whole tree
          fastBM(tree2)->y2
          fastAnc(tree2,y2)[as.character(N2)]->a
          
          ### simulating a phenotype for the binded clade by setting the same 
          ### sig2 as the positive DR case
          fastBM(Cnew.m,a=a,sig2=sigg)->ynew
          y2[tips(tree2,N2)]<-ynew
          y2->y2.list[[i]]
          
          ### computing diversifiction rates for the whole tree and comparing them to
          ### the DR of the original tree
          as.matrix(1/evol.distinct(tree2)[,2])->DRm
          rownames(DRm)<-evol.distinct(tree2)[,1]
          mean(DRt[match(Ctar.m$tip.label,rownames(DRt))])/mean(DRt[-match(Ctar.m$tip.label,rownames(DRt))])->rtm
          mean(DRm[match(Cnew.m$tip.label,rownames(DRm))])/mean(DRm[-match(Cnew.m$tip.label,rownames(DRm))])->rm
          rm/rtm->DRratio.m[i]
          
        })->test.loop
        if(class(test.loop)!="try-error") break
      })
      
      distNodes(tree1,N)[1:Nnode(tree1),1]->dn1
      distNodes(tree2,N2)[1:Nnode(tree2),1]->dn2
      length(tips(tree2,N2))-1->ll
      
      ### simulating a phenotype according to the Brownian motion for both 
      ### (positive DR and negative DR cases) trees to assess the Type I error rate
      fastBM(tree1)->yB1->yB.list[[i]]
      fastBM(tree2)->yB2->yBm.list[[i]]
      
      if(node){
        ### sampling a random node from both (positive DR and negative DR cases) trees
        ### to set as clade presumed to shift under Type I error testing framework
        subtrees(tree1)->subtr1
        sapply(subtr1[which(sapply(subtr1,Ntip)>10&sapply(subtr1,Ntip)<Ntip(tree1)*0.5)],
               function(x) getMRCA(tree1,x$tip.label))->Nsam
        sample(Nsam,1)->Nb1
        
        subtrees(tree2)->subtr2
        sapply(subtr2[which(sapply(subtr2,Ntip)>10&sapply(subtr2,Ntip)<Ntip(tree2)*0.5)],
               function(x) getMRCA(tree2,x$tip.label))->Nsam2
        sample(Nsam2,1)->Nb2
      }
      
      if(BMrates){
        if(auto){
          #### Analyses of data with positive DR shift ####
          ### performing search.shiftDR on tree and data to automatically locate possible 
          ### clade shifts in both DR and phenotypic evolutionary rates
          search.shiftDR(tree=tree1,y=y1,select.test="shift",f=(Ntip(Cnew)-1),nrep=100,foldername = tempdir())->SSdr->SSdr.list[[i]]
          SSdr[[3]][[2]]->res
          
          ### repeating the analyses on the BM phenotype to assess the Type I error rate
          search.shiftDR(tree=tree1,y=yB1,select.test="shift",f=(Ntip(Cnew)-1),nrep=100,foldername = tempdir())->SSdrB->SSdrB.list[[i]]
          SSdrB[[3]][[2]]->resB
          
          #### Analyses of data with negative DR shift ####
          ### performing search.shiftDR on tree and data to automatically locate possible 
          ### clade shifts in both DR and phenotypic evolutionary rates
          search.shiftDR(tree=tree2,y=y2,select.test="shift",f=ll,nrep=100,foldername = tempdir())->SSdr.m->SSdrm.list[[i]]
          SSdr.m[[3]][[2]]->res.m
          
          ### repeating the analyses on the BM phenotype to assess the Type I error rate
          search.shiftDR(tree=tree2,y=yB2,select.test="shift",f=ll,nrep=100,foldername = tempdir())->SSdrB.m->SSdrBm.list[[i]]
          SSdrB.m[[3]][[2]]->resB.m
          
          list(c(node.plus=N,node.minus=N2),
               brown.plus=data.frame(resB,dist.node=dn1[match(rownames(resB),names(dn1))]),
               brown.minus=data.frame(resB.m,dist.node=dn2[match(rownames(resB.m),names(dn2))]),
               test.plus=data.frame(res,dist.node=dn1[match(rownames(res),names(dn1))]),
               test.minus=data.frame(res.m,dist.node=dn2[match(rownames(res.m),names(dn2))]))->res.list[[i]]
          
          if(length(Soverfit$BMrates)>0){
            PH[i]<-phr<-sample(seq(0.1,0.25,0.1),1)
            Soverfit$BMrates$clus->cc
            Soverfit$BMrates$nsim->ss
            ### testing the robstuness of search.shiftDR results to sampling effect by performing overfitDR
            overfitDR(ssDR=SSdr,y=y1,test.type="double",nsim=ss,clus=cc,ph=phr,foldername=tempdir())->OFdr.list[[i]]
            overfitDR(ssDR=SSdrB,y=yB1,test.type="double",nsim=ss,clus=cc,ph=phr,foldername=tempdir())->OFdrB.list[[i]]
            overfitDR(ssDR=SSdr.m,y=y2,test.type="double",nsim=ss,clus=cc,ph=phr,foldername=tempdir())->OFdrm.list[[i]]
            overfitDR(ssDR=SSdrB.m,y=yB2,test.type="double",nsim=ss,clus=cc,ph=phr,foldername=tempdir())->OFdrBm.list[[i]]
          }
        }
        if(node){
          #### Analyses of data with positive DR shift ####
          ### performing search.shiftDR on tree and data by setting the clade expected to shift
          search.shiftDR(tree=tree1,y=y1,node=N,select.test="shift",foldername = tempdir())->SSdrN->SSdr.node.list[[i]]
          SSdrN[[3]][[1]]->resN
          
          ### repeating the analyses on a BM phenotype and random node to assess the Type I error rate
          search.shiftDR(tree=tree1,y=yB1,node=Nb1,select.test="shift",foldername = tempdir())->SSdrBN->SSdrB.node.list[[i]]
          SSdrBN[[3]][[1]]->resBN
          
          #### Analyses of data with negative DR shift ####
          ### performing search.shiftDR on tree and data by setting the clade expected to shift
          search.shiftDR(tree=tree2,y=y2,node=N2,select.test="shift",foldername = tempdir())->SSdrN.m->SSdrm.node.list[[i]]
          SSdrN.m[[3]][[1]]->resN.m
          
          ### repeating the analyses on a BM phenotype and random node to assess the Type I error rate
          search.shiftDR(tree=tree2,y=yB2,node=Nb2,select.test="shift",foldername = tempdir())->SSdrBmN->SSdrBm.node.list[[i]]
          SSdrBmN[[3]][[1]]->resBN.m
          
          list(brown.plus=resBN,
               brown.minus=resBN.m,
               test.plus=resN,
               test.minus=resN.m)->res.node.list[[i]]
          
          if(length(Soverfit$BMrates)>0){
            PH[i]<-phr<-sample(seq(0.1,0.25,0.1),1)
            Soverfit$BMrates$clus->cc
            Soverfit$BMrates$nsim->ss
            ### testing the robstuness of search.shiftDR results to sampling effect by performing overfitDR
            overfitDR(ssDR=SSdrN,y=y1,test.type="double",nsim=ss,clus=cc,ph=phr,foldername=tempdir())->OFdr.node.list[[i]]
            overfitDR(ssDR=SSdrBN,y=yB1,test.type="double",nsim=ss,clus=cc,ph=phr,foldername=tempdir())->OFdrB.node.list[[i]]
            overfitDR(ssDR=SSdrN.m,y=y2,test.type="double",nsim=ss,clus=cc,ph=phr,foldername=tempdir())->OFdrm.node.list[[i]]
            overfitDR(ssDR=SSdrBmN,y=yB2,test.type="double",nsim=ss,clus=cc,ph=phr,foldername=tempdir())->OFdrBm.node.list[[i]]
          }
        }
      }
      
      if(RRrates){
        ### performing RRphylo on tree and data pairs
        RRphylo(tree1,y1)->RR1 
        RRphylo(tree1,yB1)->RRb1
        RRphylo(tree2,y2)->RR2
        RRphylo(tree2,yB2)->RRb2
        
        if(auto){ 
          #### Analyses of data with positive DR shift ####
          ### performing search.shiftDR on RR object and data to automatically locate possible 
          ### clade shifts in both DR and RR rates
          search.shiftDR(RR=RR1,y=y1,select.test="shift",f=(Ntip(Cnew)-1),nrep=100,foldername = tempdir())->SSdrRR->SSdrRR.list[[i]]
          SSdrRR[[3]][[2]]->resRR
          
          ### repeating the analyses on the BM phenotype to assess the Type I error rate
          search.shiftDR(RR=RRb1,y=yB1,select.test="shift",f=(Ntip(Cnew)-1),nrep=100,foldername = tempdir())->SSdrRRB->SSdrBRR.list[[i]]
          SSdrRRB[[3]][[2]]->resBRR
          
          
          #### Analyses of data with negative DR shift ####
          ### search.shiftDR on RR object and data to automatically locate possible 
          ### clade shifts in both DR and RR rates
          search.shiftDR(RR=RR2,y=y2,select.test="shift",f=ll,nrep=100,foldername = tempdir())->SSdrRR.m->SSdrmRR.list[[i]]
          SSdrRR.m[[3]][[2]]->resRR.m
          
          ### repeating the analyses on the BM phenotype to assess the Type I error rate
          search.shiftDR(RR=RRb2,y=yB2,select.test="shift",f=ll,nrep=100,foldername = tempdir())->SSdrRRBm->SSdrBmRR.list[[i]]
          SSdrRRBm[[3]][[2]]->resBRR.m
          
          list(c(node.plus=N,node.minus=N2),
               brown.plus=data.frame(resBRR,dist.node=dn1[match(rownames(resBRR),names(dn1))]),
               brown.minus=data.frame(resBRR.m,dist.node=dn2[match(rownames(resBRR.m),names(dn2))]),
               test.plus=data.frame(resRR,dist.node=dn1[match(rownames(resRR),names(dn1))]),
               test.minus=data.frame(resRR.m,dist.node=dn2[match(rownames(resRR.m),names(dn2))]))->resRR.list[[i]]
          
          if(length(Soverfit$RRrates)>0){
            PH[i]<-phr<-sample(seq(0.1,0.25,0.1),1)
            Soverfit$BMrates$clus->cc
            Soverfit$BMrates$nsim->ss
            ### testing the robstuness of search.shiftDR results to sampling effect by performing overfitDR
            overfitDR(ssDR=SSdrRR,y=y1,RR=RR1,test.type="double",nsim=ss,clus=cc,ph=phr,foldername=tempdir())->OFdrRR.list[[i]]
            overfitDR(ssDR=SSdrRRB,y=yB1,RR=RRb1,test.type="double",nsim=ss,clus=cc,ph=phr,foldername=tempdir())->OFdrBRR.list[[i]]
            overfitDR(ssDR=SSdrRR.m,y=y2,RR=RR2,test.type="double",nsim=ss,clus=cc,ph=phr,foldername=tempdir())->OFdrmRR.list[[i]]
            overfitDR(ssDR=SSdrRRBm,y=yB2,RR=RRb2,test.type="double",nsim=ss,clus=cc,ph=phr,foldername=tempdir())->OFdrBmRR.list[[i]]
          }
        }
        if(node){
          #### Analyses of data with positive DR shift ####
          ### performing search.shiftDR on RR object and data by setting the clade expected to shift
          search.shiftDR(RR=RR1,y=y1,node=N,select.test="shift",foldername = tempdir())->SSdrRRN->SSdrRR.node.list[[i]]
          SSdrRRN[[3]][[1]]->resNRR
          
          ### repeating the analyses on BM phenotype and random node to assess the Type I error rate
          search.shiftDR(RR=RRb1,y=yB1,node=Nb1,select.test="shift",foldername = tempdir())->SSdrRRBN->SSdrBRR.node.list[[i]]
          SSdrRRBN[[3]][[1]]->resBNRR
          
          
          #### Analyses of data with negative DR shift ####
          ### performing search.shiftDR on RR object and data by setting the clade expected to shift
          search.shiftDR(RR=RR2,y=y2,node=N2,select.test="shift",foldername = tempdir())->SSdrRRN.m->SSdrmRR.node.list[[i]]
          SSdrRRN.m[[3]][[1]]->resN.mRR
          
          ### simulating a phenotype according to the Brownian motion for the whole tree
          ### and repeating the analyses to assess the Type I error rate
          search.shiftDR(RR=RRb2,y=yB2,node=Nb2,select.test="shift",foldername = tempdir())->SSdrRRBmN->SSdrBmRR.node.list[[i]]
          SSdrRRBmN[[3]][[1]]->resBN.mRR
          
          list(brown.plus=resBNRR,
               brown.minus=resBN.mRR,
               test.plus=resNRR,
               test.minus=resN.mRR)->resRR.node.list[[i]]
          
          if(length(Soverfit$RRrates)>0){
            PH[i]<-phr<-sample(seq(0.1,0.25,0.1),1)
            Soverfit$BMrates$clus->cc
            Soverfit$BMrates$nsim->ss
            ### testing the robstuness of search.shiftDR results to sampling effect by performing overfitDR
            overfitDR(ssDR=SSdrRRN,RR=RR1,y=y1,test.type="double",nsim=ss,clus=cc,ph=phr,foldername=tempdir())->OFdrRR.node.list[[i]]
            overfitDR(ssDR=SSdrRRBN,RR=RRb1,y=yB1,test.type="double",nsim=ss,clus=cc,ph=phr,foldername=tempdir())->OFdrBRR.node.list[[i]]
            overfitDR(ssDR=SSdrRRN.m,RR=RR2,y=y2,test.type="double",nsim=ss,clus=cc,ph=phr,foldername=tempdir())->OFdrmRR.node.list[[i]]
            overfitDR(ssDR=SSdrRRBmN,RR=RRb2,y=yB2,test.type="double",nsim=ss,clus=cc,ph=phr,foldername=tempdir())->OFdrBmRR.node.list[[i]]
          }
        }
      }
      
      ### plotting the trees with colored nodes pertaining to binded clades
      layout(1:2)
      plot(tree1,show.tip.label = F,no.margin = T)->tp
      legend(x=0,y=max(tp$y.lim),legend=i,bty="n",text.font=2,text.col="red",xjust=0.5)
      nodelabels(bg="w",frame="n",col="red",node=N)
      if(BMrates) if(!is.null(SSdr$shift$single.clades)) nodelabels(bg="w",frame="n",col="green",node=as.numeric(rownames(SSdr$shift$single.clades)))
      if(RRrates) if(!is.null(SSdrRR$shift$single.clades)) nodelabels(bg="w",frame="n",col="purple",node=as.numeric(rownames(SSdrRR$shift$single.clades)))
      
      plot(tree2,show.tip.label = F,no.margin = T)
      nodelabels(bg="w",frame="n",col="red",node=N2)
      if(BMrates) if(!is.null(SSdr.m$shift$single.clades)) nodelabels(bg="w",frame="n",col="green",node=as.numeric(rownames(SSdr.m$shift$single.clades)))
      if(RRrates) if(!is.null(SSdrRR.m$shift$single.clades)) nodelabels(bg="w",frame="n",col="purple",node=as.numeric(rownames(SSdrRR.m$shift$single.clades)))
      
    }
    
    ##### Character-associated diversification #####
    if(!is.null(trend.args)){
      ### creating a phylogenetic tree with > 180 species
      repeat({
        sim.bdtree(b = .5, d = 0.2)->tree
        if(Ntip(tree)>180) break
      })
      
      ### simulating a phenotype according to the Brownian motion
      fastBM(tree)->y
      
      ### sampling the parameters for the species sampling
      sample(5:10,1)->div->Div.trend[i] ### the fraction of tips to be sampled
      sample(seq(1,10,0.01),1)->esp->Esp.trend[i] ### the exponent conditioning the sampling probability of y values
      
      ### sampling the species where new tips are to be binded
      ### for i < reps/2 lower y values have higher sampling probability weights;
      ### for i > reps/2 higher y values have higher sampling probability weights;
      if(i<round(reps/2)) names(y[sample(1:length(y),Ntip(tree)/div,prob=(length(y)-rank(y))^esp)])->samtips else names(y[sample(1:length(y),Ntip(tree)/div,prob=rank(y)^esp)])->samtips
      
      tree->treeO
      y->yO
      
      ### binding a new tip halfway the leaf of each sampled tip
      for(w in 1:length(samtips)){
        tree$edge[,2]->eds
        eds[which(eds<=Ntip(tree))]<-tree$tip.label[eds[which(eds<=Ntip(tree))]]
        tree$edge.length[which(eds==samtips[w])]/2->pos
        bind.tip(tree,paste("t",w,sep=""),edge.length = pos,
                 position=pos,where=which(tree$tip.label==samtips[w]))->tree
        
      }
      
      ### reconstructing phenotypic values at the binded tips according to the Brownian motion
      fastBM(tree)->xx
      y[match(names(xx),names(y))]->yy
      names(xx)[which(is.na(match(names(xx),names(y))))]->names(yy)[which(is.na(match(names(xx),names(y))))]
      data.frame(species=names(yy),yy,xx)->yy
      phylopars(yy,tree)->res
      res$anc_recon[1:Ntip(tree),1]->y.full1
      
      
      ### repeating the procedure for the newly generated tree
      if(i<round(reps/2)) names(y.full1[sample(1:length(y.full1),Ntip(tree)/div,prob=(length(y.full1)-rank(y.full1))^esp)])->samtips else 
        names(y.full1[sample(1:length(y.full1),Ntip(tree)/div,prob=rank(y.full1)^esp)])->samtips
      
      for(w in 1:length(samtips)){
        tree$edge[,2]->eds
        eds[which(eds<=Ntip(tree))]<-tree$tip.label[eds[which(eds<=Ntip(tree))]]
        tree$edge.length[which(eds==samtips[w])]/2->pos
        bind.tip(tree,paste("x",w,sep=""),edge.length = pos,
                 position=pos,where=which(tree$tip.label==samtips[w]))->tree
      }
      Ntip(tree)/Ntip(treeO)->tree.ratio.trend[i]
      tree->tree.trend.list[[i]]
      
      ### reconstructing phenotypic values at the binded tips according to the Brownian motion
      fastBM(tree)->xx
      y.full1[match(names(xx),names(y.full1))]->yy
      names(xx)[which(is.na(match(names(xx),names(y.full1))))]->names(yy)[which(is.na(match(names(xx),names(y.full1))))]
      data.frame(species=names(yy),yy,xx)->yy
      phylopars(yy,tree)->res
      res$anc_recon[1:Ntip(tree),1]->y.full->y.trend.list[[i]]
      ecdf_fun(y,mean(y.full))->means.shift.trend[i]
      
      ### plotting the tree with colored leaves pertaining to binded tips
      colo<-rep("black",length(tree$edge[,2]))
      names(colo)<-tree$edge[,2]
      names(colo)[which(as.numeric(names(colo))<=Ntip(tree))]<-tree$tip.label
      colo[which(names(colo)%in%samtips)]<-"red"
      colo[grep("t",names(colo))]<-"green"
      colo[grep("x",names(colo))]<-"blue"
      plot(tree,show.tip.label = FALSE,no.margin=TRUE,edge.color = colo)->tp
      legend(x=0,y=max(tp$y.lim),legend=i,bty="n",text.font=2,text.col="red",xjust=0.5)
      
      
      ### performing search.shiftDR on tree and data to test the association between DR and phenotypic values
      search.shiftDR(tree,y=y.full,select.test="trend",foldername = tempdir())->SSdrT->SSdr.trend.list[[i]]
      if(trend.args$ESsim) essim(tree,y.full)->ESsims->ESsims.list[[i]]
      
      ### simulating a phenotype according to the Brownian motion for the whole tree
      ### and repeating the analysis to assess the Type I error rate
      fastBM(tree)->yB
      yB->yB.trend.list[[i]]
      search.shiftDR(tree,y=yB,select.test="trend",foldername=tempdir())->SSdrBT->SSdrB.trend.list[[i]]
      if(trend.args$ESsim) essim(tree,yB)->ESsimsB->ESsimsB.list[[i]]
      
      rbind(brown=SSdrBT$trend.res,test=SSdrT$trend.res)->res.trend.list[[i]]
      if(trend.args$ESsim) rbind(brown=ESsimsB,test=ESsims)->res.ESsim.list[[i]]
      
      if(length(trend.args$overfit)>0){
        ### testing the robustness of search.shiftDR results to sampling effects by applying overfitDR
        PH.trend[i]<-phr<-sample(seq(0.2,0.5,0.1),1)
        trend.args$overfit$clus->cc
        trend.args$overfit$nsim->ss
        overfitDR(ssDR=SSdrT,y=y.full,test.type="double",nsim=ss,clus=cc,ph=phr,foldername = tempdir())->OFdrT->OFdr.trend.list[[i]]
        overfitDR(ssDR=SSdrBT,y=yB,test.type="double",nsim=ss,clus=cc,ph=phr,foldername = tempdir())->OFdrBT->OFdrB.trend.list[[i]]
        rbind(brown=OFdrBT$trend,test=OFdrT$trend)->resOF.trend.list[[i]]
      }
    }
    
    i=i+1
  }
  
  if(!is.null(trend.args)){
    list(tree.trend.list,y.trend.list,yB.trend.list,
         res.trend.list,resOF.trend.list,
         PH.trend,Div.trend,Esp.trend,tree.ratio.trend,means.shift.trend,
         SSdr.trend.list,SSdrB.trend.list,
         OFdr.trend.list,OFdrB.trend.list,
         ESsims.list,ESsimsB.list,res.ESsim.list)->finalres
    names(finalres)<-c("tree.trend.list","y.trend.list","yB.trend.list",
                       "res.trend.list","resOF.trend.list",
                       "PH.trend","Div.trend","Esp.trend","tree.ratio.trend","means.shift.trend",
                       "SSdr.trend.list","SSdrB.trend.list",
                       "OFdr.trend.list","OFdrB.trend.list",
                       "ESsims.list","ESsimsB.list","res.ESsim.list")

    if(any(lapply(finalres,length)==0)) finalres[-(which(lapply(finalres,length)==0))]->finalres

    list2env(finalres,envir=.GlobalEnv)
    
    as.data.frame(do.call(rbind,lapply(res.trend.list[1:round(reps/2)],function(x) x[1,])))->ty1m.trend
    as.data.frame(do.call(rbind,lapply(res.trend.list[1:round(reps/2)],function(x) x[2,])))->ty2m.trend
    as.data.frame(do.call(rbind,lapply(res.trend.list[(round(reps/2)+1):reps],function(x) x[1,])))->ty1p.trend
    as.data.frame(do.call(rbind,lapply(res.trend.list[(round(reps/2)+1):reps],function(x) x[2,])))->ty2p.trend
    
    rbind(c(power=length(which(ty2p.trend$p.random>=0.95)),
            misclass=length(which(ty2p.trend$p.random<=0.05)),
            type1=length(which(ty1p.trend$p.random<=0.025|ty1p.trend$p.random>=0.975)))/nrow(ty2p.trend)*100,
          c(power=length(which(ty2m.trend$p.random<=0.05)),
            misclass=length(which(ty2m.trend$p.random>=0.95)),
            type1=length(which(ty1m.trend$p.random<=0.025|ty1m.trend$p.random>=0.975)))/nrow(ty1m.trend)*100)->trend.res
    rownames(trend.res)<-c("positive","negative")
    
    if(trend.args$ESsim){
      as.data.frame(do.call(rbind,lapply(res.ESsim.list[1:round(reps/2)],function(x) x[1,])))->ty1m.ESsim
      as.data.frame(do.call(rbind,lapply(res.ESsim.list[1:round(reps/2)],function(x) x[2,])))->ty2m.ESsim
      as.data.frame(do.call(rbind,lapply(res.ESsim.list[(round(reps/2)+1):reps],function(x) x[1,])))->ty1p.ESsim
      as.data.frame(do.call(rbind,lapply(res.ESsim.list[(round(reps/2)+1):reps],function(x) x[2,])))->ty2p.ESsim
      
      rbind(c(power=length(which(ty2p.ESsim[,1]>0&ty2p.ESsim[,2]<=0.05)),
              misclass=length(which(ty2p.ESsim[,1]<0&ty2p.ESsim[,2]<=0.05)),
              type1=length(which(ty1p.ESsim[,2]<=0.05)))/nrow(ty2p.ESsim)*100,
            c(power=length(which(ty2m.ESsim[,1]<0&ty2m.ESsim[,2]<=0.05)),
              misclass=length(which(ty2m.ESsim[,1]>0&ty2m.ESsim[,2]<=0.05)),
              type1=length(which(ty1m.ESsim[,2]<=0.05)))/nrow(ty2m.ESsim)*100)->ESsim.res
      rownames(ESsim.res)<-c("positive","negative")
      list(SSdr=trend.res,ESsim=ESsim.res)->trend.res
    }
  }else trend.res<-NULL
  
  if(!is.null(shift.args)){
    list(tree.list,tree2.list,y.list,y2.list,yB.list,yBm.list,
         NN,NN2,sigmas,DRratio,DRratio.m,PH,
         res.list,resRR.list,res.node.list,resRR.node.list,
         SSdr.list,SSdrB.list,SSdrm.list,SSdrBm.list,
         SSdrRR.list,SSdrBRR.list,SSdrmRR.list,SSdrBmRR.list,
         SSdr.node.list,SSdrB.node.list,SSdrm.node.list,SSdrBm.node.list,
         SSdrRR.node.list,SSdrBRR.node.list,SSdrmRR.node.list,SSdrBmRR.node.list,
         OFdr.list,OFdrB.list,OFdrm.list,OFdrBm.list,
         OFdrRR.list,OFdrBRR.list,OFdrmRR.list,OFdrBmRR.list,
         OFdr.node.list,OFdrB.node.list,OFdrm.node.list,OFdrBm.node.list,
         OFdrRR.node.list,OFdrBRR.node.list,OFdrmRR.node.list,OFdrBmRR.node.list)->finalres
    names(finalres)<-c("tree.list","tree2.list","y.list","y2.list","yB.list","yBm.list",
                       "NN","NN2","sigmas","DRratio","DRratio.m","PH",
                       "res.list","resRR.list","res.node.list","resRR.node.list",
                       "SSdr.list","SSdrB.list","SSdrm.list","SSdrBm.list",
                       "SSdrRR.list","SSdrBRR.list","SSdrmRR.list","SSdrBmRR.list",
                       "SSdr.node.list","SSdrB.node.list","SSdrm.node.list","SSdrBm.node.list",
                       "SSdrRR.node.list","SSdrBRR.node.list","SSdrmRR.node.list","SSdrBmRR.node.list",
                       "OFdr.list","OFdrB.list","OFdrm.list","OFdrBm.list",
                       "OFdrRR.list","OFdrBRR.list","OFdrmRR.list","OFdrBmRR.list",
                       "OFdr.node.list","OFdrB.node.list","OFdrm.node.list","OFdrBm.node.list",
                       "OFdrRR.node.list","OFdrBRR.node.list","OFdrmRR.node.list","OFdrBmRR.node.list")
    
    if(any(lapply(finalres,length)==0)) finalres[-(which(lapply(finalres,length)==0))]->finalres
    
    list2env(finalres,envir=.GlobalEnv)
    
    c(DRratio,DRratio.m)->DRratio.tot
    c(sigmas,sigmas)->sigma.tot
    
    if(BMrates){
      if(auto){
        lapply(1:length(SSdr.list),function(x) {
          res.list[[x]][[1]][1]->nx
          res.list[[x]][[1]][2]->nx2
          getDescendants(tree.list[[x]],nx)->des
          des[-which(des<=Ntip(tree.list[[x]]))]->des
          
          getDescendants(tree2.list[[x]],nx2)->des2
          des2[-which(des2<=Ntip(tree2.list[[x]]))]->des2
          
          sapply(rownames(res.list[[x]][[2]]),function(w) {
            if(w%in%c(getMommy(tree.list[[x]],nx),des,nx))
              round(length(tips(tree.list[[x]],as.numeric(w)))/length(tips(tree.list[[x]],nx)),3) else 0
          })->inx
          data.frame(iter=rep(x,nrow(res.list[[x]][[2]])),res.list[[x]][[2]],inx,
                     DRratio=rep(DRratio[x],nrow(res.list[[x]][[2]])),sigma=rep(sigmas[x],nrow(res.list[[x]][[2]])),
                     db=rep(nrow(res.list[[x]][[2]]),nrow(res.list[[x]][[2]])))->res.list[[x]][[2]]
          
          sapply(rownames(res.list[[x]][[3]]),function(w) {
            if(w%in%c(getMommy(tree2.list[[x]],nx2),des2,nx2))
              round(length(tips(tree2.list[[x]],as.numeric(w)))/length(tips(tree2.list[[x]],nx2)),3) else 0
          })->inx
          data.frame(iter=rep(x,nrow(res.list[[x]][[3]])),res.list[[x]][[3]],inx,
                     DRratio=rep(DRratio.m[x],nrow(res.list[[x]][[3]])),sigma=rep(sigmas[x],nrow(res.list[[x]][[3]])),
                     db=rep(nrow(res.list[[x]][[3]]),nrow(res.list[[x]][[3]])))->res.list[[x]][[3]]
          
          sapply(rownames(res.list[[x]][[4]]),function(w) {
            if(w%in%c(getMommy(tree.list[[x]],nx),des,nx))
              round(length(tips(tree.list[[x]],as.numeric(w)))/length(tips(tree.list[[x]],nx)),3) else 0
          })->inx
          data.frame(iter=rep(x,nrow(res.list[[x]][[4]])),res.list[[x]][[4]],inx,
                     DRratio=rep(DRratio[x],nrow(res.list[[x]][[4]])),sigma=rep(sigmas[x],nrow(res.list[[x]][[4]])),
                     db=rep(nrow(res.list[[x]][[4]]),nrow(res.list[[x]][[4]])))->res.list[[x]][[4]]
          
          sapply(rownames(res.list[[x]][[5]]),function(w) {
            if(w%in%c(getMommy(tree2.list[[x]],nx2),des2,nx2))
              round(length(tips(tree2.list[[x]],as.numeric(w)))/length(tips(tree2.list[[x]],nx2)),3) else 0
          })->inx
          data.frame(iter=rep(x,nrow(res.list[[x]][[5]])),res.list[[x]][[5]],inx,
                     DRratio=rep(DRratio.m[x],nrow(res.list[[x]][[5]])),sigma=rep(sigmas[x],nrow(res.list[[x]][[5]])),
                     db=rep(nrow(res.list[[x]][[5]]),nrow(res.list[[x]][[5]])))->res.list[[x]][[5]]
          
          
          list(node=res.list[[x]][[1]],brown.plus=res.list[[x]][[2]],brown.min=res.list[[x]][[3]],test.plus=res.list[[x]][[4]],test.min=res.list[[x]][[5]])
        })->res.list1
        
        c(lapply(res.list1, "[[",4),lapply(res.list1, "[[",5))->ty2
        lapply(ty2,function(x){
          if(nrow(x)>1&any(x$inx>0)) x[which(x$inx>0),] else if(nrow(x)>1&!any(x$inx>0)) x[1,] else if(nrow(x)==1) x
        })->ty2
       
        ty2[which(sapply(ty2,is.null))]<-list(ty2[which(!sapply(ty2,is.null))][[1]][-c(1:nrow(ty2[which(!sapply(ty2,is.null))][[1]])),])
         
        ty2[which(DRratio.tot>1&sigma.tot>=1)]->ty2.pp
        ty2[which(DRratio.tot>1&sigma.tot<1)]->ty2.pm
        ty2[which(DRratio.tot<1&sigma.tot>=1)]->ty2.mp
        ty2[which(DRratio.tot<1&sigma.tot<1)]->ty2.mm
        
        c(lapply(res.list1, "[[",2),lapply(res.list1, "[[",3))->ty1
        ty1[which(DRratio.tot>1&sigma.tot>=1)]->ty1.pp
        ty1[which(DRratio.tot>1&sigma.tot<1)]->ty1.pm
        ty1[which(DRratio.tot<1&sigma.tot>=1)]->ty1.mp
        ty1[which(DRratio.tot<1&sigma.tot<1)]->ty1.mm
        
        
        rbind(round(c(empty=length(which(sapply(ty2.pp,nrow)==0))*100/length(ty2.pp),
                      double=nrow(subset(do.call(rbind,ty2.pp),inx>0&db>1))*100/nrow(do.call(rbind,ty2.pp)),
                      error=nrow(subset(do.call(rbind,ty2.pp),inx==0))*100/nrow(do.call(rbind,ty2.pp)),
                      DRpBRp=nrow(subset(do.call(rbind,ty2.pp),inx>0&DR.diff>0&sigma.node>sigma.back))*100/nrow(do.call(rbind,ty2.pp)),
                      DRpBRm=nrow(subset(do.call(rbind,ty2.pp),inx>0&DR.diff>0&sigma.node<sigma.back))*100/nrow(do.call(rbind,ty2.pp)),
                      DRmBRp=nrow(subset(do.call(rbind,ty2.pp),inx>0&DR.diff<0&sigma.node>sigma.back))*100/nrow(do.call(rbind,ty2.pp)),
                      DRmBRm=nrow(subset(do.call(rbind,ty2.pp),inx>0&DR.diff<0&sigma.node<sigma.back))*100/nrow(do.call(rbind,ty2.pp)),
                      type1=nrow(do.call(rbind,ty1.pp))*100/length(ty1.pp)),3),
              round(c(empty=length(which(sapply(ty2.pm,nrow)==0))*100/length(ty2.pm),
                      double=nrow(subset(do.call(rbind,ty2.pm),inx>0&db>1))*100/nrow(do.call(rbind,ty2.pm)),
                      error=nrow(subset(do.call(rbind,ty2.pm),inx==0))*100/nrow(do.call(rbind,ty2.pm)),
                      DRpBRp=nrow(subset(do.call(rbind,ty2.pm),inx>0&DR.diff>0&sigma.node>sigma.back))*100/nrow(do.call(rbind,ty2.pm)),
                      DRpBRm=nrow(subset(do.call(rbind,ty2.pm),inx>0&DR.diff>0&sigma.node<sigma.back))*100/nrow(do.call(rbind,ty2.pm)),
                      DRmBRp=nrow(subset(do.call(rbind,ty2.pm),inx>0&DR.diff<0&sigma.node>sigma.back))*100/nrow(do.call(rbind,ty2.pm)),
                      DRmBRm=nrow(subset(do.call(rbind,ty2.pm),inx>0&DR.diff<0&sigma.node<sigma.back))*100/nrow(do.call(rbind,ty2.pm)),
                      type1=nrow(do.call(rbind,ty1.pm))*100/length(ty1.pm)),3),
              round(c(empty=length(which(sapply(ty2.mp,nrow)==0))*100/length(ty2.mp),
                      double=nrow(subset(do.call(rbind,ty2.mp),inx>0&db>1))*100/nrow(do.call(rbind,ty2.mp)),
                      error=nrow(subset(do.call(rbind,ty2.mp),inx==0))*100/nrow(do.call(rbind,ty2.mp)),
                      DRpBRp=nrow(subset(do.call(rbind,ty2.mp),inx>0&DR.diff>0&sigma.node>sigma.back))*100/nrow(do.call(rbind,ty2.mp)),
                      DRpBRm=nrow(subset(do.call(rbind,ty2.mp),inx>0&DR.diff>0&sigma.node<sigma.back))*100/nrow(do.call(rbind,ty2.mp)),
                      DRmBRp=nrow(subset(do.call(rbind,ty2.mp),inx>0&DR.diff<0&sigma.node>sigma.back))*100/nrow(do.call(rbind,ty2.mp)),
                      DRmBRm=nrow(subset(do.call(rbind,ty2.mp),inx>0&DR.diff<0&sigma.node<sigma.back))*100/nrow(do.call(rbind,ty2.mp)),
                      type1=nrow(do.call(rbind,ty1.mp))*100/length(ty1.mp)),3),
              round(c(empty=length(which(sapply(ty2.mm,nrow)==0))*100/length(ty2.mm),
                      double=nrow(subset(do.call(rbind,ty2.mm),inx>0&db>1))*100/nrow(do.call(rbind,ty2.mm)),
                      error=nrow(subset(do.call(rbind,ty2.mm),inx==0))*100/nrow(do.call(rbind,ty2.mm)),
                      DRpBRp=nrow(subset(do.call(rbind,ty2.mm),inx>0&DR.diff>0&sigma.node>sigma.back))*100/nrow(do.call(rbind,ty2.mm)),
                      DRpBRm=nrow(subset(do.call(rbind,ty2.mm),inx>0&DR.diff>0&sigma.node<sigma.back))*100/nrow(do.call(rbind,ty2.mm)),
                      DRmBRp=nrow(subset(do.call(rbind,ty2.mm),inx>0&DR.diff<0&sigma.node>sigma.back))*100/nrow(do.call(rbind,ty2.mm)),
                      DRmBRm=nrow(subset(do.call(rbind,ty2.mm),inx>0&DR.diff<0&sigma.node<sigma.back))*100/nrow(do.call(rbind,ty2.mm)),
                      type1=nrow(do.call(rbind,ty1.mm))*100/length(ty1.mm)),3))->restabBM
        rownames(restabBM)<-c("DRpBRp","DRpBRm","DRmBRp","DRmBRm")
      }else restabBM<-NULL
      
      if(node){
        rbind(data.frame(iter=seq(1,length(SSdr.node.list)),do.call(rbind,lapply(res.node.list,"[[",1)),DRratio=DRratio,sigma=sigmas,Blim=sapply(lapply(SSdrB.node.list,"[[",3),"[[",3),sim="pos"),
              data.frame(iter=seq(1,length(SSdr.node.list)),do.call(rbind,lapply(res.node.list,"[[",2)),DRratio=DRratio.m,sigma=sigmas,Blim=sapply(lapply(SSdrBm.node.list,"[[",3),"[[",3),sim="neg"))->ty1.node
        rbind(data.frame(iter=seq(1,length(SSdr.node.list)),do.call(rbind,lapply(res.node.list,"[[",3)),DRratio=DRratio,sigma=sigmas,Blim=sapply(lapply(SSdr.node.list,"[[",3),"[[",3),sim="pos"),
              data.frame(iter=seq(1,length(SSdr.node.list)),do.call(rbind,lapply(res.node.list,"[[",4)),DRratio=DRratio.m,sigma=sigmas,Blim=sapply(lapply(SSdrm.node.list,"[[",3),"[[",3),sim="neg"))->ty2.node
        ty2.node$tip.ratio<-1
        ty1.node$tip.ratio<-1
        ty2.node$inx<-NA
        ty1.node$inx<-NA
        
        if(any(ty2.node$p.DR>=0.975&ty2.node$p.sigma<=0.05&ty2.node$sigma.ratio>=ty2.node$Blim|
               ty2.node$p.DR<=0.025&ty2.node$p.sigma<=0.05&ty2.node$sigma.ratio>=ty2.node$Blim))
          ty2.node[which(ty2.node$p.DR>=0.975&ty2.node$p.sigma<=0.05&ty2.node$sigma.ratio>=ty2.node$Blim|
                           ty2.node$p.DR<=0.025&ty2.node$p.sigma<=0.05&ty2.node$sigma.ratio>=ty2.node$Blim),]$inx<-1
        if(any(is.na(ty2.node$inx))) ty2.node$inx[which(is.na(ty2.node$inx))]<-0
        
        if(any(ty1.node$p.DR>=0.975&ty1.node$p.sigma<=0.05&ty1.node$sigma.ratio>=ty1.node$Blim|
                 ty1.node$p.DR<=0.025&ty1.node$p.sigma<=0.05&ty1.node$sigma.ratio>=ty1.node$Blim))
          ty1.node[which(ty1.node$p.DR>=0.975&ty1.node$p.sigma<=0.05&ty1.node$sigma.ratio>=ty1.node$Blim|
                         ty1.node$p.DR<=0.025&ty1.node$p.sigma<=0.05&ty1.node$sigma.ratio>=ty1.node$Blim),]$inx<-1
        if(any(is.na(ty1.node$inx))) ty1.node$inx[which(is.na(ty1.node$inx))]<-0
        
        #### path nodes ####
        rbind(do.call(rbind,lapply(1:length(SSdr.node.list),function(w){
          SSdr.node.list[[w]][[3]]->x
          if(!is.null(x[[2]])) {
            do.call(rbind,lapply(1:nrow(x[[2]]),function(q)
              data.frame(iter=w,x[[2]][q,],DRratio=DRratio[w],sigma=sigmas[w],Blim=x[[3]],sim="pos",
                         tip.ratio=length(tips(tree.list[[w]],as.numeric(rownames(x[[2]])[q])))/length(tips(tree.list[[w]],NN[w])),
                         inx=0.5)))->aaa
            aaa[which.min(aaa$tip.ratio-1),]
          }
        })),
        do.call(rbind,lapply(1:length(SSdr.node.list),function(w){  
          SSdrm.node.list[[w]][[3]]->x
          if(!is.null(x[[2]])) {
            do.call(rbind,lapply(1:nrow(x[[2]]),function(q)
              data.frame(iter=w,x[[2]][q,],DRratio=DRratio.m[w],sigma=sigmas[w],Blim=x[[3]],sim="neg",
                         tip.ratio=length(tips(tree2.list[[w]],as.numeric(rownames(x[[2]])[q])))/length(tips(tree2.list[[w]],NN2[w])),
                         inx=0.5)))->aaa
            aaa[which.min(aaa$tip.ratio-1),]
          }
          
        })))->ty2.path
        
        if(any(ty2.node$inx==0)){
          ty2.node[which(ty2.node$inx==0),]->ty2.nosig
          for(i in 1:nrow(ty2.nosig)){
            if(ty2.nosig[i,]$sim=="pos") 
              if(ty2.nosig[i,]$iter%in%ty2.path[which(ty2.path$sim=="pos"),]$iter) ty2.path[which(ty2.path$sim=="pos"),][which(ty2.path[which(ty2.path$sim=="pos"),]$iter==ty2.nosig[i,]$iter),]->ty2.nosig[i,]
            
            if(ty2.nosig[i,]$sim=="neg") 
              if(ty2.nosig[i,]$iter%in%ty2.path[which(ty2.path$sim=="neg"),]$iter) ty2.path[which(ty2.path$sim=="neg"),][which(ty2.path[which(ty2.path$sim=="neg"),]$iter==ty2.nosig[i,]$iter),]->ty2.nosig[i,]
            
          }
          ty2.nosig->ty2.node[which(ty2.node$inx==0),]
        }
        
        rbind(do.call(rbind,lapply(1:length(SSdr.node.list),function(w){
          SSdrB.node.list[[w]][[3]]->x
          if(!is.null(x[[2]])) {
            do.call(rbind,lapply(1:nrow(x[[2]]),function(q)
              data.frame(iter=w,x[[2]][q,],DRratio=DRratio[w],sigma=sigmas[w],Blim=x[[3]],sim="pos",
                         tip.ratio=length(tips(tree.list[[w]],as.numeric(rownames(x[[2]])[q])))/length(tips(tree.list[[w]],NN[w])),
                         inx=0.5)))->aaa
            aaa[which.min(aaa$tip.ratio-1),]
          }
        })),
        do.call(rbind,lapply(1:length(SSdr.node.list),function(w){  
          SSdrBm.node.list[[w]][[3]]->x
          if(!is.null(x[[2]])) {
            do.call(rbind,lapply(1:nrow(x[[2]]),function(q)
              data.frame(iter=w,x[[2]][q,],DRratio=DRratio.m[w],sigma=sigmas[w],Blim=x[[3]],sim="neg",
                         tip.ratio=length(tips(tree2.list[[w]],as.numeric(rownames(x[[2]])[q])))/length(tips(tree2.list[[w]],NN2[w])),
                         inx=0.5)))->aaa
            aaa[which.min(aaa$tip.ratio-1),]
          }
          
        })))->ty1.path
        
        if(any(ty1.node$inx==0)){
          ty1.node[which(ty1.node$inx==0),]->ty1.nosig
          for(i in 1:nrow(ty1.nosig)){
            if(ty1.nosig[i,]$sim=="pos") 
              if(ty1.nosig[i,]$iter%in%ty1.path[which(ty1.path$sim=="pos"),]$iter) ty1.path[which(ty1.path$sim=="pos"),][which(ty1.path[which(ty1.path$sim=="pos"),]$iter==ty1.nosig[i,]$iter),]->ty1.nosig[i,]
            
            if(ty1.nosig[i,]$sim=="neg") 
              if(ty1.nosig[i,]$iter%in%ty1.path[which(ty1.path$sim=="neg"),]$iter) ty1.path[which(ty1.path$sim=="neg"),][which(ty1.path[which(ty1.path$sim=="neg"),]$iter==ty1.nosig[i,]$iter),]->ty1.nosig[i,]
            
          }
          ty1.nosig->ty1.node[which(ty1.node$inx==0),]
        }
        
        subset(ty2.node,DRratio>1&sigma>=1)->ty2.node.pp
        subset(ty2.node,DRratio>1&sigma<1)->ty2.node.pm
        subset(ty2.node,DRratio<1&sigma>=1)->ty2.node.mp
        subset(ty2.node,DRratio<1&sigma<1)->ty2.node.mm
        
        subset(ty1.node,DRratio>1&sigma>=1)->ty1.node.pp
        subset(ty1.node,DRratio>1&sigma<1)->ty1.node.pm
        subset(ty1.node,DRratio<1&sigma>=1)->ty1.node.mp
        subset(ty1.node,DRratio<1&sigma<1)->ty1.node.mm
        
        
        rbind(round(c(error=nrow(subset(ty2.node.pp,inx<1))*100/nrow(ty2.node.pp),
                      DRpBRp=nrow(subset(ty2.node.pp,inx==1&DR.diff>0&sigma.node>sigma.back))*100/nrow(ty2.node.pp),
                      DRpBRm=nrow(subset(ty2.node.pp,inx==1&DR.diff>0&sigma.node<sigma.back))*100/nrow(ty2.node.pp),
                      DRmBRp=nrow(subset(ty2.node.pp,inx==1&DR.diff<0&sigma.node>sigma.back))*100/nrow(ty2.node.pp),
                      DRmBRm=nrow(subset(ty2.node.pp,inx==1&DR.diff<0&sigma.node<sigma.back))*100/nrow(ty2.node.pp),
                      type1=length(which(ty1.node.pp$inx>0))*100/nrow(ty1.node.pp)),3),
              round(c(error=nrow(subset(ty2.node.pm,inx<1))*100/nrow(ty2.node.pm),
                      DRpBRp=nrow(subset(ty2.node.pm,inx==1&DR.diff>0&sigma.node>sigma.back))*100/nrow(ty2.node.pm),
                      DRpBRm=nrow(subset(ty2.node.pm,inx==1&DR.diff>0&sigma.node<sigma.back))*100/nrow(ty2.node.pm),
                      DRmBRp=nrow(subset(ty2.node.pm,inx==1&DR.diff<0&sigma.node>sigma.back))*100/nrow(ty2.node.pm),
                      DRmBRm=nrow(subset(ty2.node.pm,inx==1&DR.diff<0&sigma.node<sigma.back))*100/nrow(ty2.node.pm),
                      type1=length(which(ty1.node.pm$inx>0))*100/nrow(ty1.node.pm)),3),
              round(c(error=nrow(subset(ty2.node.mp,inx<1))*100/nrow(ty2.node.mp),
                      DRpBRp=nrow(subset(ty2.node.mp,inx==1&DR.diff>0&sigma.node>sigma.back))*100/nrow(ty2.node.mp),
                      DRpBRm=nrow(subset(ty2.node.mp,inx==1&DR.diff>0&sigma.node<sigma.back))*100/nrow(ty2.node.mp),
                      DRmBRp=nrow(subset(ty2.node.mp,inx==1&DR.diff<0&sigma.node>sigma.back))*100/nrow(ty2.node.mp),
                      DRmBRm=nrow(subset(ty2.node.mp,inx==1&DR.diff<0&sigma.node<sigma.back))*100/nrow(ty2.node.mp),
                      type1=length(which(ty1.node.mp$inx>0))*100/nrow(ty1.node.mp)),3),
              round(c(error=nrow(subset(ty2.node.mm,inx<1))*100/nrow(ty2.node.mm),
                      DRpBRp=nrow(subset(ty2.node.mm,inx==1&DR.diff>0&sigma.node>sigma.back))*100/nrow(ty2.node.mm),
                      DRpBRm=nrow(subset(ty2.node.mm,inx==1&DR.diff>0&sigma.node<sigma.back))*100/nrow(ty2.node.mm),
                      DRmBRp=nrow(subset(ty2.node.mm,inx==1&DR.diff<0&sigma.node>sigma.back))*100/nrow(ty2.node.mm),
                      DRmBRm=nrow(subset(ty2.node.mm,inx==1&DR.diff<0&sigma.node<sigma.back))*100/nrow(ty2.node.mm),
                      type1=length(which(ty1.node.mm$inx>0))*100/nrow(ty1.node.mm)),3))->restab.focal.nodeBM
        
        rbind(round(c(DRpBRp.path=nrow(subset(ty2.node.pp,inx==0.5&DR.diff>0&sigma.node>sigma.back))*100/nrow(ty2.node.pp),
                      DRpBRm.path=nrow(subset(ty2.node.pp,inx==0.5&DR.diff>0&sigma.node<sigma.back))*100/nrow(ty2.node.pp),
                      DRmBRp.path=nrow(subset(ty2.node.pp,inx==0.5&DR.diff<0&sigma.node>sigma.back))*100/nrow(ty2.node.pp),
                      DRmBRm.path=nrow(subset(ty2.node.pp,inx==0.5&DR.diff<0&sigma.node<sigma.back))*100/nrow(ty2.node.pp),
                      tip.ratio={if(length(subset(ty2.node.pp,inx==0.5)$tip.ratio)>0) mean(subset(ty2.node.pp,inx==0.5)$tip.ratio) else 0},
                      type1=nrow(subset(ty1.node.pp,inx==0.5))*100/nrow(ty1.node.pp)),3),
              round(c(DRpBRp.path=nrow(subset(ty2.node.pm,inx==0.5&DR.diff>0&sigma.node>sigma.back))*100/nrow(ty2.node.pm),
                      DRpBRm.path=nrow(subset(ty2.node.pm,inx==0.5&DR.diff>0&sigma.node<sigma.back))*100/nrow(ty2.node.pm),
                      DRmBRp.path=nrow(subset(ty2.node.pm,inx==0.5&DR.diff<0&sigma.node>sigma.back))*100/nrow(ty2.node.pm),
                      DRmBRm.path=nrow(subset(ty2.node.pm,inx==0.5&DR.diff<0&sigma.node<sigma.back))*100/nrow(ty2.node.pm),
                      tip.ratio={if(length(subset(ty2.node.pm,inx==0.5)$tip.ratio)>0) mean(subset(ty2.node.pm,inx==0.5)$tip.ratio) else 0},
                      type1=nrow(subset(ty1.node.pm,inx==0.5))*100/nrow(ty1.node.pm)),3),
              round(c(DRpBRp.path=nrow(subset(ty2.node.mp,inx==0.5&DR.diff>0&sigma.node>sigma.back))*100/nrow(ty2.node.mp),
                      DRpBRm.path=nrow(subset(ty2.node.mp,inx==0.5&DR.diff>0&sigma.node<sigma.back))*100/nrow(ty2.node.mp),
                      DRmBRp.path=nrow(subset(ty2.node.mp,inx==0.5&DR.diff<0&sigma.node>sigma.back))*100/nrow(ty2.node.mp),
                      DRmBRm.path=nrow(subset(ty2.node.mp,inx==0.5&DR.diff<0&sigma.node<sigma.back))*100/nrow(ty2.node.mp),
                      tip.ratio={if(length(subset(ty2.node.mp,inx==0.5)$tip.ratio)>0) mean(subset(ty2.node.mp,inx==0.5)$tip.ratio) else 0},
                      type1=nrow(subset(ty1.node.mp,inx==0.5))*100/nrow(ty1.node.mp)),3),
              round(c(DRpBRp.path=nrow(subset(ty2.node.mm,inx==0.5&DR.diff>0&sigma.node>sigma.back))*100/nrow(ty2.node.mm),
                      DRpBRm.path=nrow(subset(ty2.node.mm,inx==0.5&DR.diff>0&sigma.node<sigma.back))*100/nrow(ty2.node.mm),
                      DRmBRp.path=nrow(subset(ty2.node.mm,inx==0.5&DR.diff<0&sigma.node>sigma.back))*100/nrow(ty2.node.mm),
                      DRmBRm.path=nrow(subset(ty2.node.mm,inx==0.5&DR.diff<0&sigma.node<sigma.back))*100/nrow(ty2.node.mm),
                      tip.ratio={if(length(subset(ty2.node.mm,inx==0.5)$tip.ratio)>0) mean(subset(ty2.node.mm,inx==0.5)$tip.ratio) else 0},
                      type1=nrow(subset(ty1.node.mm,inx==0.5))*100/nrow(ty1.node.mm)),3))->restab.path.nodeBM
        rownames(restab.focal.nodeBM)<-rownames(restab.path.nodeBM)<-c("DRpBRp","DRpBRm","DRmBRp","DRmBRm")
      }else {
        restab.focal.nodeBM<-NULL
        restab.path.nodeBM<-NULL
      }
      
      list(automatic=restabBM,node=list(focal=restab.focal.nodeBM,path=restab.path.nodeBM))->shiftBMres
      
    }else shiftBMres<-NULL
    
    if(RRrates){
      if(auto){
        lapply(1:length(SSdrRR.list),function(x) {
          resRR.list[[x]][[1]][1]->nx
          resRR.list[[x]][[1]][2]->nx2
          getDescendants(tree.list[[x]],nx)->des
          des[-which(des<=Ntip(tree.list[[x]]))]->des
          
          getDescendants(tree2.list[[x]],nx2)->des2
          des2[-which(des2<=Ntip(tree2.list[[x]]))]->des2
          
          sapply(rownames(resRR.list[[x]][[2]]),function(w) {
            if(w%in%c(getMommy(tree.list[[x]],nx),des,nx))
              round(length(tips(tree.list[[x]],as.numeric(w)))/length(tips(tree.list[[x]],nx)),3) else 0
          })->inx
          data.frame(iter=rep(x,nrow(resRR.list[[x]][[2]])),resRR.list[[x]][[2]],inx,
                     DRratio=rep(DRratio[x],nrow(resRR.list[[x]][[2]])),sigma=rep(sigmas[x],nrow(resRR.list[[x]][[2]])),
                     db=rep(nrow(resRR.list[[x]][[2]]),nrow(resRR.list[[x]][[2]])))->resRR.list[[x]][[2]]
          
          sapply(rownames(resRR.list[[x]][[3]]),function(w) {
            if(w%in%c(getMommy(tree2.list[[x]],nx2),des2,nx2))
              round(length(tips(tree2.list[[x]],as.numeric(w)))/length(tips(tree2.list[[x]],nx2)),3) else 0
          })->inx
          data.frame(iter=rep(x,nrow(resRR.list[[x]][[3]])),resRR.list[[x]][[3]],inx,
                     DRratio=rep(DRratio.m[x],nrow(resRR.list[[x]][[3]])),sigma=rep(sigmas[x],nrow(resRR.list[[x]][[3]])),
                     db=rep(nrow(resRR.list[[x]][[3]]),nrow(resRR.list[[x]][[3]])))->resRR.list[[x]][[3]]
          
          sapply(rownames(resRR.list[[x]][[4]]),function(w) {
            if(w%in%c(getMommy(tree.list[[x]],nx),des,nx))
              round(length(tips(tree.list[[x]],as.numeric(w)))/length(tips(tree.list[[x]],nx)),3) else 0
          })->inx
          data.frame(iter=rep(x,nrow(resRR.list[[x]][[4]])),resRR.list[[x]][[4]],inx,
                     DRratio=rep(DRratio[x],nrow(resRR.list[[x]][[4]])),sigma=rep(sigmas[x],nrow(resRR.list[[x]][[4]])),
                     db=rep(nrow(resRR.list[[x]][[4]]),nrow(resRR.list[[x]][[4]])))->resRR.list[[x]][[4]]
          
          sapply(rownames(resRR.list[[x]][[5]]),function(w) {
            if(w%in%c(getMommy(tree2.list[[x]],nx2),des2,nx2))
              round(length(tips(tree2.list[[x]],as.numeric(w)))/length(tips(tree2.list[[x]],nx2)),3) else 0
          })->inx
          data.frame(iter=rep(x,nrow(resRR.list[[x]][[5]])),resRR.list[[x]][[5]],inx,
                     DRratio=rep(DRratio.m[x],nrow(resRR.list[[x]][[5]])),sigma=rep(sigmas[x],nrow(resRR.list[[x]][[5]])),
                     db=rep(nrow(resRR.list[[x]][[5]]),nrow(resRR.list[[x]][[5]])))->resRR.list[[x]][[5]]
          
          
          list(node=resRR.list[[x]][[1]],brown.plus=resRR.list[[x]][[2]],brown.min=resRR.list[[x]][[3]],test.plus=resRR.list[[x]][[4]],test.min=resRR.list[[x]][[5]])
        })->res.list1RR
        
        c(lapply(res.list1RR, "[[",4),lapply(res.list1RR, "[[",5))->ty2RR
        
        lapply(ty2RR,function(x){
          if(nrow(x)>1&any(x$inx>0)) x[which(x$inx>0),] else if(nrow(x)>1&!any(x$inx>0)) x[1,] else if(nrow(x)==1) x
        })->ty2RR
        
        ty2RR[which(sapply(ty2RR,is.null))]<-list(ty2RR[which(!sapply(ty2RR,is.null))][[1]][-c(1:nrow(ty2RR[which(!sapply(ty2RR,is.null))][[1]])),])
        
        ty2RR[which(DRratio.tot>1&sigma.tot>=1)]->ty2RR.pp
        ty2RR[which(DRratio.tot>1&sigma.tot<1)]->ty2RR.pm
        ty2RR[which(DRratio.tot<1&sigma.tot>=1)]->ty2RR.mp
        ty2RR[which(DRratio.tot<1&sigma.tot<1)]->ty2RR.mm
        
        c(lapply(res.list1RR, "[[",2),lapply(res.list1RR, "[[",3))->ty1RR
        ty1RR[which(DRratio.tot>1&sigma.tot>=1)]->ty1RR.pp
        ty1RR[which(DRratio.tot>1&sigma.tot<1)]->ty1RR.pm
        ty1RR[which(DRratio.tot<1&sigma.tot>=1)]->ty1RR.mp
        ty1RR[which(DRratio.tot<1&sigma.tot<1)]->ty1RR.mm
        
        
        rbind(round(c(empty=length(which(sapply(ty2RR.pp,nrow)==0))*100/length(ty2RR.pp),
                      double=nrow(subset(do.call(rbind,ty2RR.pp),inx>0&db>1))*100/nrow(do.call(rbind,ty2RR.pp)),
                      error=nrow(subset(do.call(rbind,ty2RR.pp),inx==0))*100/nrow(do.call(rbind,ty2RR.pp)),
                      DRpRRp=nrow(subset(do.call(rbind,ty2RR.pp),inx>0&DR.diff>0&phen.diff>0))*100/nrow(do.call(rbind,ty2RR.pp)),
                      DRpRRm=nrow(subset(do.call(rbind,ty2RR.pp),inx>0&DR.diff>0&phen.diff<0))*100/nrow(do.call(rbind,ty2RR.pp)),
                      DRmRRp=nrow(subset(do.call(rbind,ty2RR.pp),inx>0&DR.diff<0&phen.diff>0))*100/nrow(do.call(rbind,ty2RR.pp)),
                      DRmRRm=nrow(subset(do.call(rbind,ty2RR.pp),inx>0&DR.diff<0&phen.diff<0))*100/nrow(do.call(rbind,ty2RR.pp)),
                      type1=nrow(do.call(rbind,ty1RR.pp))*100/length(ty1RR.pp)),3),
              round(c(empty=length(which(sapply(ty2RR.pm,nrow)==0))*100/length(ty2RR.pm),
                      double=nrow(subset(do.call(rbind,ty2RR.pm),inx>0&db>1))*100/nrow(do.call(rbind,ty2RR.pm)),
                      error=nrow(subset(do.call(rbind,ty2RR.pm),inx==0))*100/nrow(do.call(rbind,ty2RR.pm)),
                      DRpRRp=nrow(subset(do.call(rbind,ty2RR.pm),inx>0&DR.diff>0&phen.diff>0))*100/nrow(do.call(rbind,ty2RR.pm)),
                      DRpRRm=nrow(subset(do.call(rbind,ty2RR.pm),inx>0&DR.diff>0&phen.diff<0))*100/nrow(do.call(rbind,ty2RR.pm)),
                      DRmRRp=nrow(subset(do.call(rbind,ty2RR.pm),inx>0&DR.diff<0&phen.diff>0))*100/nrow(do.call(rbind,ty2RR.pm)),
                      DRmRRm=nrow(subset(do.call(rbind,ty2RR.pm),inx>0&DR.diff<0&phen.diff<0))*100/nrow(do.call(rbind,ty2RR.pm)),
                      type1=nrow(do.call(rbind,ty1RR.pm))*100/length(ty1RR.pm)),3),
              round(c(empty=length(which(sapply(ty2RR.mp,nrow)==0))*100/length(ty2RR.mp),
                      double=nrow(subset(do.call(rbind,ty2RR.mp),inx>0&db>1))*100/nrow(do.call(rbind,ty2RR.mp)),
                      error=nrow(subset(do.call(rbind,ty2RR.mp),inx==0))*100/nrow(do.call(rbind,ty2RR.mp)),
                      DRpRRp=nrow(subset(do.call(rbind,ty2RR.mp),inx>0&DR.diff>0&phen.diff>0))*100/nrow(do.call(rbind,ty2RR.mp)),
                      DRpRRm=nrow(subset(do.call(rbind,ty2RR.mp),inx>0&DR.diff>0&phen.diff<0))*100/nrow(do.call(rbind,ty2RR.mp)),
                      DRmRRp=nrow(subset(do.call(rbind,ty2RR.mp),inx>0&DR.diff<0&phen.diff>0))*100/nrow(do.call(rbind,ty2RR.mp)),
                      DRmRRm=nrow(subset(do.call(rbind,ty2RR.mp),inx>0&DR.diff<0&phen.diff<0))*100/nrow(do.call(rbind,ty2RR.mp)),
                      type1=nrow(do.call(rbind,ty1RR.mp))*100/length(ty1RR.mp)),3),
              round(c(empty=length(which(sapply(ty2RR.mm,nrow)==0))*100/length(ty2RR.mm),
                      double=nrow(subset(do.call(rbind,ty2RR.mm),inx>0&db>1))*100/nrow(do.call(rbind,ty2RR.mm)),
                      error=nrow(subset(do.call(rbind,ty2RR.mm),inx==0))*100/nrow(do.call(rbind,ty2RR.mm)),
                      DRpRRp=nrow(subset(do.call(rbind,ty2RR.mm),inx>0&DR.diff>0&phen.diff>0))*100/nrow(do.call(rbind,ty2RR.mm)),
                      DRpRRm=nrow(subset(do.call(rbind,ty2RR.mm),inx>0&DR.diff>0&phen.diff<0))*100/nrow(do.call(rbind,ty2RR.mm)),
                      DRmRRp=nrow(subset(do.call(rbind,ty2RR.mm),inx>0&DR.diff<0&phen.diff>0))*100/nrow(do.call(rbind,ty2RR.mm)),
                      DRmRRm=nrow(subset(do.call(rbind,ty2RR.mm),inx>0&DR.diff<0&phen.diff<0))*100/nrow(do.call(rbind,ty2RR.mm)),
                      type1=nrow(do.call(rbind,ty1RR.mm))*100/length(ty1RR.mm)),3))->restabRR
        rownames(restabRR)<-c("DRpRRp","DRpRRm","DRmRRp","DRmRRm")
      }else restabRR<-NULL
      
      if(node){
        rbind(data.frame(iter=seq(1,length(SSdrRR.node.list)),do.call(rbind,lapply(resRR.node.list,"[[",1)),DRratio=DRratio,sigma=sigmas,Blim=sapply(lapply(SSdrBRR.node.list,"[[",3),"[[",3),sim="pos"),
              data.frame(iter=seq(1,length(SSdrRR.node.list)),do.call(rbind,lapply(resRR.node.list,"[[",2)),DRratio=DRratio.m,sigma=sigmas,Blim=sapply(lapply(SSdrBmRR.node.list,"[[",3),"[[",3),sim="neg"))->ty1RR.node
        rbind(data.frame(iter=seq(1,length(SSdrRR.node.list)),do.call(rbind,lapply(resRR.node.list,"[[",3)),DRratio=DRratio,sigma=sigmas,Blim=sapply(lapply(SSdrRR.node.list,"[[",3),"[[",3),sim="pos"),
              data.frame(iter=seq(1,length(SSdrRR.node.list)),do.call(rbind,lapply(resRR.node.list,"[[",4)),DRratio=DRratio.m,sigma=sigmas,Blim=sapply(lapply(SSdrmRR.node.list,"[[",3),"[[",3),sim="neg"))->ty2RR.node
        ty2RR.node$tip.ratio<-1
        ty1RR.node$tip.ratio<-1
        ty2RR.node$inx<-NA
        ty1RR.node$inx<-NA
        
        if(any(ty2RR.node$p.DR>=0.975&ty2RR.node$p.phen<=0.025&ty2RR.node$sigma.ratio>=ty2RR.node$Blim|
               ty2RR.node$p.DR<=0.025&ty2RR.node$p.phen<=0.025&ty2RR.node$sigma.ratio>=ty2RR.node$Blim|
               ty2RR.node$p.DR>=0.975&ty2RR.node$p.phen>=0.975&ty2RR.node$sigma.ratio>=ty2RR.node$Blim|
               ty2RR.node$p.DR<=0.025&ty2RR.node$p.phen>=0.975&ty2RR.node$sigma.ratio>=ty2RR.node$Blim))
          ty2RR.node[which(ty2RR.node$p.DR>=0.975&ty2RR.node$p.phen<=0.025&ty2RR.node$sigma.ratio>=ty2RR.node$Blim|
                             ty2RR.node$p.DR<=0.025&ty2RR.node$p.phen<=0.025&ty2RR.node$sigma.ratio>=ty2RR.node$Blim|
                             ty2RR.node$p.DR>=0.975&ty2RR.node$p.phen>=0.975&ty2RR.node$sigma.ratio>=ty2RR.node$Blim|
                             ty2RR.node$p.DR<=0.025&ty2RR.node$p.phen>=0.975&ty2RR.node$sigma.ratio>=ty2RR.node$Blim),]$inx<-1
        if(any(is.na(ty2RR.node$inx))) ty2RR.node$inx[which(is.na(ty2RR.node$inx))]<-0
        
        if(any(ty1RR.node$p.DR>=0.975&ty1RR.node$p.phen<=0.025&ty1RR.node$sigma.ratio>=ty1RR.node$Blim|
               ty1RR.node$p.DR<=0.025&ty1RR.node$p.phen<=0.025&ty1RR.node$sigma.ratio>=ty1RR.node$Blim|
               ty1RR.node$p.DR>=0.975&ty1RR.node$p.phen>=0.975&ty1RR.node$sigma.ratio>=ty1RR.node$Blim|
               ty1RR.node$p.DR<=0.025&ty1RR.node$p.phen>=0.975&ty1RR.node$sigma.ratio>=ty1RR.node$Blim))
          ty1RR.node[which(ty1RR.node$p.DR>=0.975&ty1RR.node$p.phen<=0.025&ty1RR.node$sigma.ratio>=ty1RR.node$Blim|
                             ty1RR.node$p.DR<=0.025&ty1RR.node$p.phen<=0.025&ty1RR.node$sigma.ratio>=ty1RR.node$Blim|
                             ty1RR.node$p.DR>=0.975&ty1RR.node$p.phen>=0.975&ty1RR.node$sigma.ratio>=ty1RR.node$Blim|
                             ty1RR.node$p.DR<=0.025&ty1RR.node$p.phen>=0.975&ty1RR.node$sigma.ratio>=ty1RR.node$Blim),]$inx<-1
        if(any(is.na(ty1RR.node$inx))) ty1RR.node$inx[which(is.na(ty1RR.node$inx))]<-0
        
        #### path nodes ####
        rbind(do.call(rbind,lapply(1:length(SSdrRR.node.list),function(w){
          SSdrRR.node.list[[w]][[3]]->x
          if(!is.null(x[[2]])) {
            do.call(rbind,lapply(1:nrow(x[[2]]),function(q)
              data.frame(iter=w,x[[2]][q,],DRratio=DRratio[w],sigma=sigmas[w],Blim=x[[3]],sim="pos",
                         tip.ratio=length(tips(tree.list[[w]],as.numeric(rownames(x[[2]])[q])))/length(tips(tree.list[[w]],NN[w])),
                         inx=0.5)))->aaa
            aaa[which.min(aaa$tip.ratio-1),]
          }
        })),
        do.call(rbind,lapply(1:length(SSdrmRR.node.list),function(w){  
          SSdrmRR.node.list[[w]][[3]]->x
          if(!is.null(x[[2]])) {
            do.call(rbind,lapply(1:nrow(x[[2]]),function(q)
              data.frame(iter=w,x[[2]][q,],DRratio=DRratio.m[w],sigma=sigmas[w],Blim=x[[3]],sim="neg",
                         tip.ratio=length(tips(tree2.list[[w]],as.numeric(rownames(x[[2]])[q])))/length(tips(tree2.list[[w]],NN2[w])),
                         inx=0.5)))->aaa
            aaa[which.min(aaa$tip.ratio-1),]
          }
          
        })))->ty2RR.path
        
        if(any(ty2RR.node$inx==0)){
          ty2RR.node[which(ty2RR.node$inx==0),]->ty2RR.nosig
          for(i in 1:nrow(ty2RR.nosig)){
            if(ty2RR.nosig[i,]$sim=="pos") 
              if(ty2RR.nosig[i,]$iter%in%ty2RR.path[which(ty2RR.path$sim=="pos"),]$iter) ty2RR.path[which(ty2RR.path$sim=="pos"),][which(ty2RR.path[which(ty2RR.path$sim=="pos"),]$iter==ty2RR.nosig[i,]$iter),]->ty2RR.nosig[i,]
            
            if(ty2RR.nosig[i,]$sim=="neg") 
              if(ty2RR.nosig[i,]$iter%in%ty2RR.path[which(ty2RR.path$sim=="neg"),]$iter) ty2RR.path[which(ty2RR.path$sim=="neg"),][which(ty2RR.path[which(ty2RR.path$sim=="neg"),]$iter==ty2RR.nosig[i,]$iter),]->ty2RR.nosig[i,]
            
          }
          ty2RR.nosig->ty2RR.node[which(ty2RR.node$inx==0),]
        }
        
        rbind(do.call(rbind,lapply(1:length(SSdrBRR.node.list),function(w){
          SSdrBRR.node.list[[w]][[3]]->x
          if(!is.null(x[[2]])) {
            do.call(rbind,lapply(1:nrow(x[[2]]),function(q)
              data.frame(iter=w,x[[2]][q,],DRratio=DRratio[w],sigma=sigmas[w],Blim=x[[3]],sim="pos",
                         tip.ratio=length(tips(tree.list[[w]],as.numeric(rownames(x[[2]])[q])))/length(tips(tree.list[[w]],NN[w])),
                         inx=0.5)))->aaa
            aaa[which.min(aaa$tip.ratio-1),]
          }
        })),
        do.call(rbind,lapply(1:length(SSdrBmRR.node.list),function(w){  
          SSdrBmRR.node.list[[w]][[3]]->x
          if(!is.null(x[[2]])) {
            do.call(rbind,lapply(1:nrow(x[[2]]),function(q)
              data.frame(iter=w,x[[2]][q,],DRratio=DRratio.m[w],sigma=sigmas[w],Blim=x[[3]],sim="neg",
                         tip.ratio=length(tips(tree2.list[[w]],as.numeric(rownames(x[[2]])[q])))/length(tips(tree2.list[[w]],NN2[w])),
                         inx=0.5)))->aaa
            aaa[which.min(aaa$tip.ratio-1),]
          }
          
        })))->ty1RR.path
        
        if(any(ty1RR.node$inx==0)){
          ty1RR.node[which(ty1RR.node$inx==0),]->ty1RR.nosig
          for(i in 1:nrow(ty1RR.nosig)){
            if(ty1RR.nosig[i,]$sim=="pos") 
              if(ty1RR.nosig[i,]$iter%in%ty1RR.path[which(ty1RR.path$sim=="pos"),]$iter) ty1RR.path[which(ty1RR.path$sim=="pos"),][which(ty1RR.path[which(ty1RR.path$sim=="pos"),]$iter==ty1RR.nosig[i,]$iter),]->ty1RR.nosig[i,]
            
            if(ty1RR.nosig[i,]$sim=="neg") 
              if(ty1RR.nosig[i,]$iter%in%ty1RR.path[which(ty1RR.path$sim=="neg"),]$iter) ty1RR.path[which(ty1RR.path$sim=="neg"),][which(ty1RR.path[which(ty1RR.path$sim=="neg"),]$iter==ty1RR.nosig[i,]$iter),]->ty1RR.nosig[i,]
            
          }
          ty1RR.nosig->ty1RR.node[which(ty1RR.node$inx==0),]
        }
        
        subset(ty2RR.node,DRratio>1&sigma>=1)->ty2RR.node.pp
        subset(ty2RR.node,DRratio>1&sigma<1)->ty2RR.node.pm
        subset(ty2RR.node,DRratio<1&sigma>=1)->ty2RR.node.mp
        subset(ty2RR.node,DRratio<1&sigma<1)->ty2RR.node.mm
        
        subset(ty1RR.node,DRratio>1&sigma>=1)->ty1RR.node.pp
        subset(ty1RR.node,DRratio>1&sigma<1)->ty1RR.node.pm
        subset(ty1RR.node,DRratio<1&sigma>=1)->ty1RR.node.mp
        subset(ty1RR.node,DRratio<1&sigma<1)->ty1RR.node.mm
        
        
        rbind(round(c(error=nrow(subset(ty2RR.node.pp,inx<1))*100/nrow(ty2RR.node.pp),
                      DRpRRp=nrow(subset(ty2RR.node.pp,inx==1&DR.diff>0&phen.diff>0))*100/nrow(ty2RR.node.pp),
                      DRpRRm=nrow(subset(ty2RR.node.pp,inx==1&DR.diff>0&phen.diff<0))*100/nrow(ty2RR.node.pp),
                      DRmRRp=nrow(subset(ty2RR.node.pp,inx==1&DR.diff<0&phen.diff>0))*100/nrow(ty2RR.node.pp),
                      DRmRRm=nrow(subset(ty2RR.node.pp,inx==1&DR.diff<0&phen.diff<0))*100/nrow(ty2RR.node.pp),
                      type1=length(which(ty1RR.node.pp$inx>0))*100/nrow(ty1RR.node.pp)),3),
              round(c(error=nrow(subset(ty2RR.node.pm,inx<1))*100/nrow(ty2RR.node.pm),
                      DRpRRp=nrow(subset(ty2RR.node.pm,inx==1&DR.diff>0&phen.diff>0))*100/nrow(ty2RR.node.pm),
                      DRpRRm=nrow(subset(ty2RR.node.pm,inx==1&DR.diff>0&phen.diff<0))*100/nrow(ty2RR.node.pm),
                      DRmRRp=nrow(subset(ty2RR.node.pm,inx==1&DR.diff<0&phen.diff>0))*100/nrow(ty2RR.node.pm),
                      DRmRRm=nrow(subset(ty2RR.node.pm,inx==1&DR.diff<0&phen.diff<0))*100/nrow(ty2RR.node.pm),
                      type1=length(which(ty1RR.node.pm$inx>0))*100/nrow(ty1RR.node.pm)),3),
              round(c(error=nrow(subset(ty2RR.node.mp,inx<1))*100/nrow(ty2RR.node.mp),
                      DRpRRp=nrow(subset(ty2RR.node.mp,inx==1&DR.diff>0&phen.diff>0))*100/nrow(ty2RR.node.mp),
                      DRpRRm=nrow(subset(ty2RR.node.mp,inx==1&DR.diff>0&phen.diff<0))*100/nrow(ty2RR.node.mp),
                      DRmRRp=nrow(subset(ty2RR.node.mp,inx==1&DR.diff<0&phen.diff>0))*100/nrow(ty2RR.node.mp),
                      DRmRRm=nrow(subset(ty2RR.node.mp,inx==1&DR.diff<0&phen.diff<0))*100/nrow(ty2RR.node.mp),
                      type1=length(which(ty1RR.node.mp$inx>0))*100/nrow(ty1RR.node.mp)),3),
              round(c(error=nrow(subset(ty2RR.node.mm,inx<1))*100/nrow(ty2RR.node.mm),
                      DRpRRp=nrow(subset(ty2RR.node.mm,inx==1&DR.diff>0&phen.diff>0))*100/nrow(ty2RR.node.mm),
                      DRpRRm=nrow(subset(ty2RR.node.mm,inx==1&DR.diff>0&phen.diff<0))*100/nrow(ty2RR.node.mm),
                      DRmRRp=nrow(subset(ty2RR.node.mm,inx==1&DR.diff<0&phen.diff>0))*100/nrow(ty2RR.node.mm),
                      DRmRRm=nrow(subset(ty2RR.node.mm,inx==1&DR.diff<0&phen.diff<0))*100/nrow(ty2RR.node.mm),
                      type1=length(which(ty1RR.node.mm$inx>0))*100/nrow(ty1RR.node.mm)),3))->restab.focal.nodeRR
        
        rbind(round(c(DRpRRp=nrow(subset(ty2RR.node.pp,inx==0.5&DR.diff>0&phen.diff>0))*100/nrow(ty2RR.node.pp),
                      DRpRRm=nrow(subset(ty2RR.node.pp,inx==0.5&DR.diff>0&phen.diff<0))*100/nrow(ty2RR.node.pp),
                      DRmRRp=nrow(subset(ty2RR.node.pp,inx==0.5&DR.diff<0&phen.diff>0))*100/nrow(ty2RR.node.pp),
                      DRmRRm=nrow(subset(ty2RR.node.pp,inx==0.5&DR.diff<0&phen.diff<0))*100/nrow(ty2RR.node.pp),
                      tip.ratio={if(length(subset(ty2RR.node.pp,inx==0.5)$tip.ratio)>0) mean(subset(ty2RR.node.pp,inx==0.5)$tip.ratio) else 0},
                      type1=nrow(subset(ty1RR.node.pp,inx==0.5))*100/nrow(ty1RR.node.pp)),3),
              round(c(DRpRRp=nrow(subset(ty2RR.node.pm,inx==0.5&DR.diff>0&phen.diff>0))*100/nrow(ty2RR.node.pm),
                      DRpRRm=nrow(subset(ty2RR.node.pm,inx==0.5&DR.diff>0&phen.diff<0))*100/nrow(ty2RR.node.pm),
                      DRmRRp=nrow(subset(ty2RR.node.pm,inx==0.5&DR.diff<0&phen.diff>0))*100/nrow(ty2RR.node.pm),
                      DRmRRm=nrow(subset(ty2RR.node.pm,inx==0.5&DR.diff<0&phen.diff<0))*100/nrow(ty2RR.node.pm),
                      tip.ratio={if(length(subset(ty2RR.node.pm,inx==0.5)$tip.ratio)>0) mean(subset(ty2RR.node.pm,inx==0.5)$tip.ratio) else 0},
                      type1=nrow(subset(ty1RR.node.pm,inx==0.5))*100/nrow(ty1RR.node.pm)),3),
              round(c(DRpRRp=nrow(subset(ty2RR.node.mp,inx==0.5&DR.diff>0&phen.diff>0))*100/nrow(ty2RR.node.mp),
                      DRpRRm=nrow(subset(ty2RR.node.mp,inx==0.5&DR.diff>0&phen.diff<0))*100/nrow(ty2RR.node.mp),
                      DRmRRp=nrow(subset(ty2RR.node.mp,inx==0.5&DR.diff<0&phen.diff>0))*100/nrow(ty2RR.node.mp),
                      DRmRRm=nrow(subset(ty2RR.node.mp,inx==0.5&DR.diff<0&phen.diff<0))*100/nrow(ty2RR.node.mp),
                      tip.ratio={if(length(subset(ty2RR.node.mp,inx==0.5)$tip.ratio)>0) mean(subset(ty2RR.node.mp,inx==0.5)$tip.ratio) else 0},
                      type1=nrow(subset(ty1RR.node.mp,inx==0.5))*100/nrow(ty1RR.node.mp)),3),
              round(c(DRpRRp=nrow(subset(ty2RR.node.mm,inx==0.5&DR.diff>0&phen.diff>0))*100/nrow(ty2RR.node.mm),
                      DRpRRm=nrow(subset(ty2RR.node.mm,inx==0.5&DR.diff>0&phen.diff<0))*100/nrow(ty2RR.node.mm),
                      DRmRRp=nrow(subset(ty2RR.node.mm,inx==0.5&DR.diff<0&phen.diff>0))*100/nrow(ty2RR.node.mm),
                      DRmRRm=nrow(subset(ty2RR.node.mm,inx==0.5&DR.diff<0&phen.diff<0))*100/nrow(ty2RR.node.mm),
                      tip.ratio={if(length(subset(ty2RR.node.mm,inx==0.5)$tip.ratio)>0) mean(subset(ty2RR.node.mm,inx==0.5)$tip.ratio) else 0},
                      type1=nrow(subset(ty1RR.node.mm,inx==0.5))*100/nrow(ty1RR.node.mm)),3))->restab.path.nodeRR
        
        rownames(restab.focal.nodeRR)<-rownames(restab.path.nodeRR)<-c("DRpRRp","DRpRRm","DRmRRp","DRmRRm")
        
      }else {
        restab.focal.nodeRR<-NULL
        restab.path.nodeRR<-NULL
      }
      list(automatic=restabRR,node=list(focal=restab.focal.nodeRR,path=restab.path.nodeRR))->shiftRRres
    }else shiftRRres<-NULL
    
    list(BMrates=shiftBMres,RRrates=shiftRRres)->shift.res
    
  }else shift.res<-NULL
  
  # if(!is.null(trend.args)){
  #   list(tree.trend.list,y.trend.list,yB.trend.list,
  #        res.trend.list,resOF.trend.list,
  #        PH.trend,Div.trend,Esp.trend,tree.ratio.trend,means.shift.trend,
  #        SSdr.trend.list,SSdrB.trend.list,
  #        OFdr.trend.list,OFdrB.trend.list)->finalres
  #   names(finalres)<-c("tree.trend.list","y.trend.list","yB.trend.list",
  #                      "res.trend.list","resOF.trend.list",
  #                      "PH.trend","Div.trend","Esp.trend","tree.ratio.trend","means.shift.trend",
  #                      "SSdr.trend.list","SSdrB.trend.list",
  #                      "OFdr.trend.list","OFdrB.trend.list")
  # 
  #   if(any(lapply(finalres,length)==0)) finalres[-(which(lapply(finalres,length)==0))]->finalres
  # 
  #   list2env(finalres,envir=.GlobalEnv)
  # 
  #   as.data.frame(do.call(rbind,lapply(res.trend.list[1:round(reps/2)],function(x) x[1,])))->ty1m.trend
  #   as.data.frame(do.call(rbind,lapply(res.trend.list[1:round(reps/2)],function(x) x[2,])))->ty2m.trend
  #   as.data.frame(do.call(rbind,lapply(res.trend.list[(round(reps/2)+1):reps],function(x) x[1,])))->ty1p.trend
  #   as.data.frame(do.call(rbind,lapply(res.trend.list[(round(reps/2)+1):reps],function(x) x[2,])))->ty2p.trend
  # 
  #   rbind(c(power=length(which(ty2p.trend$p.random>=0.95)),
  #           misclass=length(which(ty2p.trend$p.random<=0.05)),
  #           type1=length(which(ty1p.trend$p.random<=0.025|ty1p.trend$p.random>=0.975)))/nrow(ty2p.trend)*100,
  #         c(power=length(which(ty2m.trend$p.random<=0.05)),
  #           misclass=length(which(ty2m.trend$p.random>=0.95)),
  #           type1=length(which(ty1m.trend$p.random<=0.025|ty1m.trend$p.random>=0.975)))/nrow(ty1p.trend)*100)->trend.res
  #   rownames(trend.res)<-c("positive","negative")
  # 
  # }else trend.res<-NULL
  
  return(list(summary.shift.res=shift.res,summary.trend.res=trend.res))
}

