#' Testing for clade-associated and character-associated diversification 

#' usage
#' search.shiftDR(tree=NULL,RR=NULL,y,select.test="both",node = NULL,
#' nrep = 1000, f = NULL,DR.log=NULL,foldername)
#' 
#' description 
#' The function search.shiftDR recognizes instances of clade-associated and/or character-associated diversification occurring on a phylogeny. Clade-associated diversification is meant as significant shift in taxonomic and phenotypic diversfication rates occurring on the same clade. Character-associated diversification is meant as linear relationship between diversification rates and phenotypic values.   
#' tree: a phylogenetic tree. It must be supplied when the rate of phenotypic evolution are to be computed by applying compare.evol.rates from the package geomorph (see details). The tree needs not to be ultrametric or fully dichotomous, yet polytomies may deeply influence DR values. We suggest to resolve polytomies by using the function fix.poly.
#' RR: an object generated by the function RRphylo. It must be supplied when the rates of phenotypic evolution are to be computed by applying RRphylo (see details).
#' y: either a single vector variable or a multivariate dataset. In any case, y must be a named vector.
#' select.test: whether to test for clade-associated diversification (select.test = "shift"), character-associated diversification (select.test = "trend"), or both (select.test = "both"). Is is set as "both" by default.  
#' node: the node (clade) to be tested for the rate shift. When multiple nodes are tested, they need to be written as in the example below. If 'node' is left unspecified, the function performs under the 'automatic' feature, meaning it will automatically test individual clades for deviation of their rates from the background rate of the rest of the tree (see details).
#' nrep: the number of simulations to be performed for the rate shift test, by default 'nrep' is set at 1000.
#' f: the size of the smallest clade to be tested. By default, nodes subtending to one tenth of the tree tips are tested.
#' DR.log: logical specifying whether diversification rate (DR) should be log-transformed (DR.log = TRUE) or used as they are (DR.log = FALSE) when testing for character-associated diversification. If 'DR.log' is left unspecified, the normality of DR values is assessed by the Shapiro-Wilk test, and possibily the log-transformation is applied.
#' foldername: the path of the folder where plots are to be stored.

#' details 
#' The function returns a list object including:
#' $tree the phylogenetic tree
#' $DR.log a logical specifying whether DR values have been log-transformed .
#' $shift.res the result list for clade-associated diversification testing
#' $all.clades results of shift testing for all selected nodes, regardless of significance. The dataframe always includes the difference between the mean DR of the clade and the mean DR of the rest of the tree (DR.diff), the significance level for DR.diff (p.DR), and the maximum to minimum ratio between the brownian rates computed for the clade and the rest of the tree (sigma.ratio). If rates of phenotypic evolution were computed by package geomorph's compare.evol.rates, the dataframe also includes the brownian rate of the clade (sgima.node) and the rest of the tree (sigma.back), and the significance level for the sigma.ratio (p.sigma). If RRphylo rates were used as rates of phenotypic evolution, the dataframe includes the difference between the mean absolute RR rates of the clade and the same figure for the rest of the tree (phen.diff), and the significance level for phen.diff (p.phen).
#' $single.clades results of shift testing restricted to significant clades having the largest/smallest mean DR values along a single path (see details). The dataframe is the same as described for '$all.clades'.
#' $others results of shift testing restricted to clades showing significant p-values for both DR and phenotypic evolutionary rates shifts, but having the sigma.ratio smaller than the 'max.s2.ratio.bm'.
#' $max.s2ratio.bm the maximum brownian rate value among the selected clades computed by compare.evol.rates. It is used as a threshold to select clades showing significant rate shift (see details).
#' $trend.res results of character-associated diversification. This include the slope and significance level of DR versus phenotype regression (slope and p.real respectively), and the p-value obtained by comparing the real slope to a family of 'nrep' slopes derived by simulating the phenotypes according to the BM (p.random).
#' 



search.shiftDR<-function(tree=NULL,
                         RR=NULL,
                         y,
                         select.test="both",
                         node=NULL,
                         nrep=1000,
                         f=NULL,
                         DR.log=NULL,
                         foldername){
  
  
  require(phytools)
  require(geiger)
  require(scales)
  require(geomorph)
  require(RRphylo)
  require(picante)
  require(plotrix)
  require(RColorBrewer)
  
  if(!is.null(tree)&!is.null(RR)) stop("Please indicate 'tree' to perform the analysis by using BM rate of phenotypic evolution 
                                       or 'RR' to use RRphylo rates of phenotypic evolution")
  
  if(is.null(RR)){
    if(!is.binary(tree)) warning("Tree includes polytomies, which may severely affect DR computation.
                               Consider resolving by using function fix.poly, no overfit testing will be otherwise possible")
  }else{
    tree <- RR$tree
    rates <- RR$rates
    betas<-RR$multiple.rates
  }
  
  if (is.null(nrow(y))) y <- treedata(tree, y, sort = TRUE)[[2]][, 1] else y <- treedata(tree, y, sort = TRUE)[[2]]
  
  as.matrix(1/evol.distinct(tree)[,2])->DR
  rownames(DR)<-evol.distinct(tree)[,1]
  
  #### Shift ####
  if(select.test=="both" | select.test=="shift"){
    if (is.null(node)) { 
      ##### Automatic #####
      DR->DRrates
      if(is.null(f)) f<-round(Ntip(tree)/10)
      ST <- subtrees(tree)
      len <- array()
      for (i in 1:length(ST)) len[i] <- Ntip(ST[[i]])
      st <- ST[which(len < (Ntip(tree)/2) & len > round(f))]
      node <- sapply(st, function(x) getMRCA(tree, x$tip.label))
      names(st) <- node
      
      #### Shift DR ####
      leaf2N.diff <- array()
      p.single <- array()
      for (j in 1:length(node)) {
        Ctips <- tips(tree, as.numeric(node[j]))
        leaf.rates <- DRrates[match(Ctips, rownames(DRrates)),]
        NCrates <- DRrates[-match(names(leaf.rates), rownames(DRrates)),]
        leafR <- mean(abs(leaf.rates))
        NCR <- mean(abs(NCrates))
        leaf2N.diff[j] <- leafR - NCR
        NC <- length(DRrates) - length(leaf.rates)
        C <- length(leaf.rates)
        ran.diffM <- array()
        for (i in 1:nrep) {
          ran.diffM[i] <- mean(sample(abs(DRrates), C)) - mean(sample(abs(DRrates),
                                                                      NC))
        }
        p.single[j] <- rank(c(leaf2N.diff[j], ran.diffM[1:(nrep -
                                                             1)]))[1]/nrep
      }
      names(leaf2N.diff) <- node
      names(p.single) <- node
      
      p.single->p.all### vedi i res finali
      
      if(!is.null(RR)){ #### RR rates ####
        p.single->pDR.single
        leaf2N.diff->leaf2N.DR.diff
        
        leaf2N.diff <- array()
        p.single <- array()
        for (j in 1:length(node)) {
          Cleaf <- getDescendants(tree, node[j])
          Cleaf[which(Cleaf<=Ntip(tree))]<-tree$tip.label[Cleaf[which(Cleaf<=Ntip(tree))]]
          leaf.rates <- rates[match(Cleaf, rownames(rates)),]
          NCrates <- rates[-match(names(leaf.rates), rownames(rates))]
          leaf2N.diff[j] <- mean(abs(leaf.rates)) - mean(abs(NCrates))
          NC <- length(NCrates)
          C <- length(leaf.rates)
          ran.diffM <- array()
          for (i in 1:nrep) {
            ran.diffM[i] <- mean(sample(abs(rates), C)) - mean(sample(abs(rates),
                                                                      NC))
          }
          p.single[j] <- rank(c(leaf2N.diff[j], ran.diffM[1:(nrep -
                                                               1)]))[1]/nrep
        }
        names(leaf2N.diff) <- node
        names(p.single) <- node
        
        data.frame("DR.diff"=leaf2N.DR.diff,"p.DR"=pDR.single,"phen.diff"=leaf2N.diff,"p.phen"=p.single)->allres
        
        if(length(y)>Ntip(tree)) fastBM(tree,nsim=ncol(y))->yB else fastBM(tree)->yB
        sigrat<-sigratB<-array()
        for(k in 1:nrow(allres)){
          if(length(y)>Ntip(tree)) {
            cat<-rep("a",nrow(y))
            names(cat)<-rownames(y) 
          }else{
            cat<-rep("a",length(y))
            names(cat)<-names(y)
          }
          cat[tips(tree,as.numeric(rownames(allres)[k]))]<-"nod"
          compare.evol.rates(cbind(y,y),tree,cat,print.progress = FALSE)->com
          compare.evol.rates(cbind(yB,yB),tree,cat,print.progress = FALSE)->comB
          com[[1]]->sigrat[k]
          comB[[1]]->sigratB[k]
        }
        data.frame(allres,"sigma.ratio"=sigrat)->allres
        max(sigratB)*2->Bsig
        
        allres[which(allres[,2]<=0.025|allres[,2]>=0.975),]->rescut
        if(nrow(rescut)>0) rescut[which(rescut[,4]<=0.025|rescut[,4]>=0.975),]->signres else rescut->signres
        if(nrow(signres)>0) signres[which(signres[,5]>=Bsig),]->single else signres->single
        
        if(nrow(single)>1){
          single[order(single[,3]),]->single
          ups.all<-ups <- single[which(single[,2] >= 0.975),2]
          names(ups)<-rownames(single[which(single[,2] >= 0.975),])
          dws.all<-dws <- single[which(single[,2] <= 0.025),2]
          names(dws)<-rownames(single[which(single[,2] <= 0.025),])
          
          allres.diff <- single[,3]
          names(allres.diff)<-rownames(single)
          
          if (length(dws)!=0) {
            s = 1
            repeat{
              d <- which(names(dws) %in% getDescendants(tree, names(dws)[s]))
              if (length(d) > 0) {
                allres.diff[c(match(names(dws[d]),names(allres.diff)),match(names(dws[s]),names(allres.diff)))]->cla
                names(which.max(abs(cla)))->IN
                dws[-match(names(cla[which(names(cla)!=IN)]),names(dws))]->dws
                s=1
              } else {
                dws <- dws
                s = s+1
              }
              if (s > length(dws))  break
            }
          }else dws<-NULL
          
          if (length(ups)!=0) {
            z = 1
            repeat{
              d <- which(names(ups) %in% getDescendants(tree, names(ups)[z]))
              if (length(d) > 0) {
                allres.diff[c(match(names(ups[d]),names(allres.diff)),match(names(ups[z]),names(allres.diff)))]->cla
                names(which.max(abs(cla)))->IN
                ups[-match(names(cla[which(names(cla)!=IN)]),names(ups))]->ups
                z=1
              }else{
                ups <- ups
                z = z+1
              }
              if (z > length(ups))  break
            }
          } else ups<-NULL 
          
          if(length(ups)>=1&length(dws)>=1){
            
            p.single.dws<-array()
            for (j in 1:length(dws)) {
              unlist(lapply(as.numeric(names(ups)),function(x) tips(tree,x)))->rem.ups
              Ctips.dws <- tips(tree, as.numeric(names(dws)[j]))
              leaf.rates.dws <- DRrates[match(Ctips.dws, rownames(DRrates)),]
              
              c(getMommy(tree,names(dws)[j])[1],getDescendants(tree,names(dws)[j])[1])->imm
              if(any(imm%in%names(dws.all))) imm[which(imm%in%names(dws.all))]->bar else bar<-NULL
              if(length(bar)>0){
                roll<-list(); for(w in 1:length(bar)) tips(tree,as.numeric(bar[w]))->roll[[w]]
                unique(c(unlist(roll),names(leaf.rates.dws)))->tm
                NCrates.dws <- DRrates[-c(match(tm, rownames(DRrates)),
                                          match(rem.ups, rownames(DRrates))),]
              }else{
                NCrates.dws <- DRrates[-c(match(names(leaf.rates.dws), rownames(DRrates)),
                                          match(rem.ups, rownames(DRrates))),]
              }
              
              leaf2N.diff.dws <- mean(abs(leaf.rates.dws)) - mean(abs(NCrates.dws))
              NC.dws <- length(NCrates.dws)
              C.dws <- length(leaf.rates.dws)
              
              ran.diffM.dws <- array()
              for (i in 1:nrep) {
                ran.diffM.dws[i] <- mean(sample(abs(DRrates), C.dws)) - mean(sample(abs(DRrates),NC.dws))
              }
              p.single.dws[j] <- rank(c(leaf2N.diff.dws, ran.diffM.dws[1:(nrep -1)]))[1]/nrep
            }
            
            p.single.ups <- array()
            for (j in 1:length(ups)) {
              unlist(lapply(as.numeric(names(dws)),function(x) tips(tree,x)))->rem.dws
              Ctips.ups <- tips(tree, as.numeric(names(ups)[j]))
              leaf.rates.ups <- DRrates[match(Ctips.ups, rownames(DRrates)),]
              
              c(getMommy(tree,names(ups)[j])[1],getDescendants(tree,names(ups)[j])[1])->imm
              if(any(imm%in%names(ups.all))) imm[which(imm%in%names(ups.all))]->bar else bar<-NULL
              if(length(bar)>0){
                roll<-list(); for(w in 1:length(bar)) tips(tree,as.numeric(bar[w]))->roll[[w]]
                unique(c(unlist(roll),names(leaf.rates.ups)))->tm
                NCrates.ups <- DRrates[-c(match(tm, rownames(DRrates)),
                                          match(rem.dws, rownames(DRrates))),]
              }else{
                NCrates.ups <- DRrates[-c(match(names(leaf.rates.ups), rownames(DRrates)),
                                          match(rem.dws, rownames(DRrates))),]
              }
              
              
              leaf2N.diff.ups <- mean(abs(leaf.rates.ups)) - mean(abs(NCrates.ups))
              NC.ups <- length(NCrates.ups)
              C.ups <- length(leaf.rates.ups)
              
              ran.diffM.ups <- array()
              for (i in 1:nrep) {
                ran.diffM.ups[i] <- mean(sample(abs(DRrates), C.ups)) - mean(sample(abs(DRrates),NC.ups))
              }
              p.single.ups[j] <- rank(c(leaf2N.diff.ups, ran.diffM.ups[1:(nrep -1)]))[1]/nrep
            }
            
            if(!all(c(p.single.dws>0.05,p.single.ups<0.95))){
              if(any(p.single.dws>0.05)) dws[-which(p.single.dws>0.05)]->dws
              if(any(p.single.ups<0.95)) ups[-which(p.single.ups<0.95)]->ups
            }
          }
          
          c(dws,ups)->dwup
          p.single.dwup<-array()
          for (j in 1:length(dwup)) {
            unlist(lapply(as.numeric(names(dwup)[-j]),function(x) getDescendants(tree,x)))->rem
            rem[which(rem<=Ntip(tree))]<-tree$tip.label[rem[which(rem<=Ntip(tree))]]
            Ctips.foc <- getDescendants(tree, as.numeric(names(dwup)[j]))
            Ctips.foc[which(Ctips.foc<=Ntip(tree))]<-tree$tip.label[Ctips.foc[which(Ctips.foc<=Ntip(tree))]]
            leaf.rates.foc <- rates[match(Ctips.foc, rownames(rates)),]
            NCrates.foc <- rates[-c(match(names(leaf.rates.foc), rownames(rates)),
                                    match(rem, rownames(rates))),]
            
            leaf2N.diff.foc <- mean(abs(leaf.rates.foc))-mean(abs(NCrates.foc))
            NC.foc <- length(NCrates.foc)
            C.foc <- length(leaf.rates.foc)
            
            ran.diffM.foc <- array()
            for (i in 1:nrep) {
              ran.diffM.foc[i] <- mean(sample(abs(rates), C.foc)) - mean(sample(abs(rates),NC.foc))
            }
            p.single.dwup[j] <- rank(c(leaf2N.diff.foc, ran.diffM.foc[1:(nrep -1)]))[1]/nrep
          }
          names(p.single.dwup)<-names(dwup)
          
          if(any(p.single.dwup>=0.975|p.single.dwup<=0.025)) 
            single<-single[which(rownames(single)%in%names(which(p.single.dwup>=0.975|p.single.dwup<=0.025))),] else single<-NULL
          
        }
        
        if(!is.null(single)){
          if(nrow(single)>0) signres[-match(rownames(single),rownames(signres)),]->signres else single<-NULL
        }
        
        if(nrow(signres)<1) signres<-NULL
        shift.res<-list(allres,single,signres,Bsig)
        names(shift.res)<-c("all.clades","single.clades","others","max.s2ratio.bm")
        
        pdf(file=paste(foldername, "Clade associated diversification plot.pdf",sep="/"))
        if(!is.null(shift.res$single.clades)){
          if(Ntip(tree)>100) plot(tree,edge.col="gray50", show.tip.label = FALSE) else plot(tree,edge.col="gray50", cex=.8)
          plotinfo<-get("last_plot.phylo",envir =ape::.PlotPhyloEnv)
          t(sapply(as.numeric(rownames(shift.res$single.clades)), function(x) c(plotinfo$xx[x],plotinfo$yy[x])))->xy
          rownames(xy)<-rownames(shift.res$single.clades)
          cbind(xy,NA)->xy
          
          brewer.pal(8,"Paired")[c(2,4,6,8)]->colo
          if(length(which(shift.res$single.clades$p.DR>=0.975&shift.res$single.clades$p.phen>=0.975))>0)
            xy[which(shift.res$single.clades$p.DR>=0.975&shift.res$single.clades$p.phen>=0.975),3]<-colo[3]
          if(length(which(shift.res$single.clades$p.DR>=0.975&shift.res$single.clades$p.phen<=0.025))>0)
            xy[which(shift.res$single.clades$p.DR>=0.975&shift.res$single.clades$p.phen<=0.025),3]<-colo[4]
          if(length(which(shift.res$single.clades$p.DR<=0.025&shift.res$single.clades$p.phen>=0.975))>0)
            xy[which(shift.res$single.clades$p.DR<=0.025&shift.res$single.clades$p.phen>=0.975),3]<-colo[2]
          if(length(which(shift.res$single.clades$p.DR<=0.025&shift.res$single.clades$p.phen<=0.025))>0)
            xy[which(shift.res$single.clades$p.DR<=0.025&shift.res$single.clades$p.phen<=0.025),3]<-colo[1]
          
          apply(xy,1,function(x) symbols(x[1],x[2],circles=0.2,inches=0.15,bg = alpha(x[3], 0.5), fg = x[3],add=TRUE))
          
          nodelabels(node = as.numeric(rownames(shift.res$single.clades)), adj = c(1.5,1), text = as.numeric(rownames(shift.res$single.clades)), frame = "none", bg = "white",
                     col = "black",font=2)
          
          leg<-legend("topleft",legend=c("shift in DR- and RR rate-","shift in DR- and RR rate+",
                                         "shift in DR+ and RR rate-","shift in DR+ and RR rate+"),pch=16,
                      col=colo[c(1,2,4,3)],bty="n",pt.cex=1.7,x.intersp=0.7,cex=.8,plot=FALSE)
          -0.8*leg$rect$h/(max(plotinfo$yy)-min(plotinfo$yy))->ii
          legend("topleft",legend=c("shift in DR- and RR rate-","shift in DR- and RR rate+",
                                    "shift in DR+ and RR rate-","shift in DR+ and RR rate+"),pch=16,
                 col=colo[c(1,2,4,3)],bty="n",pt.cex=1.7,x.intersp=0.7,cex=.8,inset=c(0,ii),xpd=TRUE)
        }else{
          if(Ntip(tree)>100) plot(tree, show.tip.label = FALSE) else plot(tree, cex=.8)
        }
        dev.off()
      }else{ #### BM rates ####
        if(length(y)>Ntip(tree)) fastBM(tree,nsim=ncol(y))->yB else fastBM(tree)->yB
        aovres<-list()
        sigratB<-array()
        for(k in 1:length(p.all)){
          if(length(y)>Ntip(tree)) {
            cat<-rep("a",nrow(y))
            names(cat)<-rownames(y) 
          }else{
            cat<-rep("a",length(y))
            names(cat)<-names(y)
          }
          cat[tips(tree,as.numeric(names(p.all[k])))]<-"nod"
          compare.evol.rates(cbind(y,y),tree,cat,print.progress = FALSE)->com
          c(com[[5]][which(names(com[[5]])=="nod")],com[[5]][-which(names(com[[5]])=="nod")])->sigma.g
          c(unlist(com[1:2]),sigma.g)->aovres[[k]]
          compare.evol.rates(cbind(yB,yB),tree,cat,print.progress = FALSE)[[1]]->sigratB[k]
        }
        do.call(rbind,aovres)->aovres
        rownames(aovres)<-names(p.all)
        colnames(aovres)<-c("sigma.ratio","p.sigma","sigma.node","sigma.back")
        max(sigratB)*2->Bsig
        
        data.frame("DR.diff"=leaf2N.diff,"p.DR"=p.all,aovres)->allres
        
        if(length(which(allres[,4]<=0.05))==0)  p.single<-NULL else p.single[which(allres[,4]<=0.05&allres[,3]>Bsig)]->p.single
        
        if (length(p.single[p.single >= 0.975 | p.single <=0.025])==0)  p.single<-NULL
        
        
        if (length(p.single[p.single >= 0.975 | p.single <=
                            0.025]) == 1) {
          p.single <- p.single[p.single >= 0.975 | p.single <=
                                 0.025]
          leaf2N.diff <- leaf2N.diff[match(names(p.single),
                                           names(leaf2N.diff))]
          #p.single<-NULL
        }
        
        if (length(p.single[p.single >= 0.975 | p.single <=
                            0.025]) >= 2)  {
          
          p.single <- p.single[p.single >= 0.975 | p.single <= 0.025]
          allres.diff <- allres[match(names(p.single), rownames(allres)),3]
          names(allres.diff)<-names(p.single)
          
          ups <- p.single[p.single >= 0.975]
          dws <- p.single[p.single <= 0.025]
          
          ups.all<-ups <- ups[na.omit(match(names(allres.diff[order(allres.diff,
                                                                    decreasing = FALSE)]), names(ups)))]
          dws.all<-dws <- dws[na.omit(match(names(allres.diff[order(allres.diff,
                                                                    decreasing = FALSE)]), names(dws)))]
          
          if (!is.na(mean(dws))){
            s = 1
            repeat{
              d <- which(names(dws) %in% getDescendants(tree, names(dws)[s]))
              if (length(d) > 0) {
                allres.diff[c(match(names(dws[d]),names(allres.diff)),match(names(dws[s]),names(allres.diff)))]->cla
                names(which.max(abs(cla)))->IN
                dws[-match(names(cla[which(names(cla)!=IN)]),names(dws))]->dws
                s=1
              }else{
                dws <- dws
                s = s+1
              }
              if (s > length(dws))  break
            }
          } else dws<-NULL
          
          if (!is.na(mean(ups))){
            z = 1
            repeat{
              d <- which(names(ups) %in% getDescendants(tree, names(ups)[z]))
              if (length(d) > 0) {
                allres.diff[c(match(names(ups[d]),names(allres.diff)),match(names(ups[z]),names(allres.diff)))]->cla
                names(which.max(abs(cla)))->IN
                ups[-match(names(cla[which(names(cla)!=IN)]),names(ups))]->ups
                z=1
              }else{
                ups <- ups
                z = z+1
              }
              if (z > length(ups))  break 
            }
          } else ups<-NULL
          
          if(length(ups)>=1&length(dws)>=1){
            p.single.dws<-array()
            for (j in 1:length(dws)) {
              unlist(lapply(as.numeric(names(ups)),function(x) tips(tree,x)))->rem.ups
              Ctips.dws <- tips(tree, as.numeric(names(dws)[j]))
              leaf.rates.dws <- DRrates[match(Ctips.dws, rownames(DRrates)),]
              
              c(getMommy(tree,names(dws)[j])[1],getDescendants(tree,names(dws)[j])[1])->imm
              if(any(imm%in%names(dws.all))) imm[which(imm%in%names(dws.all))]->bar else bar<-NULL
              if(length(bar)>0){
                roll<-list(); for(w in 1:length(bar)) tips(tree,as.numeric(bar[w]))->roll[[w]]
                unique(c(unlist(roll),names(leaf.rates.dws)))->tm
                NCrates.dws <- DRrates[-c(match(tm, rownames(DRrates)),
                                          match(rem.ups, rownames(DRrates))),]
              }else{
                NCrates.dws <- DRrates[-c(match(names(leaf.rates.dws), rownames(DRrates)),
                                          match(rem.ups, rownames(DRrates))),]
              }
              
              leafR.dws <- mean(abs(leaf.rates.dws))
              NCR.dws <- mean(abs(NCrates.dws))
              leaf2N.diff.dws <- leafR.dws - NCR.dws
              NC.dws <- length(NCrates.dws)
              C.dws <- length(leaf.rates.dws)
              
              ran.diffM.dws <- array()
              for (i in 1:nrep) {
                ran.diffM.dws[i] <- mean(sample(abs(DRrates), C.dws)) - mean(sample(abs(DRrates),NC.dws))
              }
              p.single.dws[j] <- rank(c(leaf2N.diff.dws, ran.diffM.dws[1:(nrep -1)]))[1]/nrep
            }
            
            p.single.ups <- array()
            for (j in 1:length(ups)) {
              unlist(lapply(as.numeric(names(dws)),function(x) tips(tree,x)))->rem.dws
              Ctips.ups <- tips(tree, as.numeric(names(ups)[j]))
              leaf.rates.ups <- DRrates[match(Ctips.ups, rownames(DRrates)),]
              
              c(getMommy(tree,names(ups)[j])[1],getDescendants(tree,names(ups)[j])[1])->imm
              if(any(imm%in%names(ups.all))) imm[which(imm%in%names(ups.all))]->bar else bar<-NULL
              if(length(bar)>0){
                roll<-list(); for(w in 1:length(bar)) tips(tree,as.numeric(bar[w]))->roll[[w]]
                unique(c(unlist(roll),names(leaf.rates.ups)))->tm
                NCrates.ups <- DRrates[-c(match(tm, rownames(DRrates)),
                                          match(rem.dws, rownames(DRrates))),]
              }else{
                NCrates.ups <- DRrates[-c(match(names(leaf.rates.ups), rownames(DRrates)),
                                          match(rem.dws, rownames(DRrates))),]
              }
              
              
              leafR.ups <- mean(abs(leaf.rates.ups))
              NCR.ups <- mean(abs(NCrates.ups))
              leaf2N.diff.ups <- leafR.ups - NCR.ups
              NC.ups <- length(NCrates.ups)
              C.ups <- length(leaf.rates.ups)
              
              ran.diffM.ups <- array()
              for (i in 1:nrep) {
                ran.diffM.ups[i] <- mean(sample(abs(DRrates), C.ups)) - mean(sample(abs(DRrates),NC.ups))
              }
              p.single.ups[j] <- rank(c(leaf2N.diff.ups, ran.diffM.ups[1:(nrep -1)]))[1]/nrep
            }
            
            
            if(!all(c(p.single.dws>0.05,p.single.ups<0.95))){
              if(any(p.single.dws>0.05)) dws[-which(p.single.dws>0.05)]->dws
              if(any(p.single.ups<0.95)) ups[-which(p.single.ups<0.95)]->ups
            }
            
          }
          
          p.single <- p.single[which(names(p.single) %in%
                                       names(c(ups, dws)))]
          
          lapply(names(p.single),function(x) tips(tree,as.numeric(x)))->sta
          names(sta)<-names(p.single)
          stack(sta)->sta1
          
          if(length(y)>Ntip(tree)){
            cat<-rep("back",nrow(y))
            names(cat)<-rownames(y)
          }else{
            cat<-rep("back",length(y))
            names(cat)<-names(y)
          }
          cat[match(sta1[,1],names(cat))]<-paste("g",sta1[,2],sep="")
          
          compare.evol.rates(cbind(y,y),tree,cat,print.progress = FALSE)->cer
          if(inherits(cer$pairwise.pvalue,"dist")) as.matrix(cer$pairwise.pvalue)[-1,1]->pcer else {
            cer$pairwise.pvalue->pcer
            names(pcer)<-cer$groups[2]
          }
          if(length(which(pcer<=0.05))>0)
            p.single <- p.single[match(gsub("g","",names(which(pcer<=0.05))),names(p.single))] else p.single<-NULL
          
          leaf2N.diff <- leaf2N.diff[match(names(p.single),
                                           names(leaf2N.diff))]
          
        }
        
        if(length(p.single)>1){
          data.frame("DR.diff"=leaf2N.diff[match(names(p.single),names(leaf2N.diff))],
                     "p.DR"=p.single,aovres[match(names(p.single),rownames(aovres)),c(1,3,4,2)])->single.clades
        }else{
          if(is.null(p.single)) single.clades<-NULL else {
            as.data.frame(matrix(c(unname(leaf2N.diff[match(names(p.single),names(leaf2N.diff))]),
                                   unname(p.single),aovres[match(names(p.single),rownames(aovres)),c(1,3,4,2)]),nrow=1))->single.clades
            
            rownames(single.clades)<-names(p.single)
            colnames(single.clades)<-c("DR.diff","p.DR",colnames(aovres)[c(1,3,4,2)])
          }
        }
        
        
        allres[,2]->resord
        if(any(resord>0.5)) 1-resord[which(resord>0.5)]->resord[which(resord>0.5)]
        allres[order(resord,allres[,4]),]->allres
        
        if(is.null(single.clades)) allres->signres else allres[-which(rownames(allres)%in%rownames(single.clades)),]->signres
        if(any(signres[,2]<=0.025&signres[,4]<=0.05|signres[,2]>=0.975&signres[,4]<=0.05))
          signres[which(signres[,2]<=0.025&signres[,4]<=0.05|signres[,2]>=0.975&signres[,4]<=0.05),]->signres else NULL->signres
        
        allres[,c(1:3,5:6,4)]->allres
        signres[,c(1:3,5:6,4)]->signres
        shift.res<-list(allres,single.clades,signres,Bsig)
        names(shift.res)<-c("all.clades","single.clades","others","max.s2ratio.bm")
        
        
        pdf(file=paste(foldername, "Clade associated diversification plot.pdf",sep="/"))
        if(!is.null(shift.res$single.clades)){
          if(Ntip(tree)>100) plot(tree,edge.col="gray50", show.tip.label = FALSE) else plot(tree,edge.col="gray50", cex=.8)
          plotinfo<-get("last_plot.phylo",envir =ape::.PlotPhyloEnv)
          t(sapply(as.numeric(rownames(shift.res$single.clades)), function(x) c(plotinfo$xx[x],plotinfo$yy[x])))->xy
          rownames(xy)<-rownames(shift.res$single.clades)
          cbind(xy,NA)->xy
          
          brewer.pal(8,"Paired")[c(2,4,6,8)]->colo
          if(length(which(shift.res$single.clades$p.DR>0.975&shift.res$single.clades$sigma.node>shift.res$single.clades$sigma.back))>0)
            xy[which(shift.res$single.clades$p.DR>0.975&shift.res$single.clades$sigma.node>shift.res$single.clades$sigma.back),3]<-colo[3]
          if(length(which(shift.res$single.clades$p.DR>0.975&shift.res$single.clades$sigma.node<shift.res$single.clades$sigma.back))>0)
            xy[which(shift.res$single.clades$p.DR>0.975&shift.res$single.clades$sigma.node<shift.res$single.clades$sigma.back),3]<-colo[4]
          if(length(which(shift.res$single.clades$p.DR<0.025&shift.res$single.clades$sigma.node>shift.res$single.clades$sigma.back))>0)
            xy[which(shift.res$single.clades$p.DR<0.025&shift.res$single.clades$sigma.node>shift.res$single.clades$sigma.back),3]<-colo[2]
          if(length(which(shift.res$single.clades$p.DR<0.025&shift.res$single.clades$sigma.node<shift.res$single.clades$sigma.back))>0)
            xy[which(shift.res$single.clades$p.DR<0.025&shift.res$single.clades$sigma.node<shift.res$single.clades$sigma.back),3]<-colo[1]
          
          apply(xy,1,function(x) symbols(x[1],x[2],circles=0.2,inches=0.15,bg = alpha(x[3], 0.5), fg = x[3],add=TRUE))
          
          nodelabels(node = as.numeric(rownames(shift.res$single.clades)), adj = c(1.5,1), text = as.numeric(rownames(shift.res$single.clades)), frame = "none", bg = "white",
                     col = "black",font=2)
          
          leg<-legend("topleft",legend=c("shift in DR- and BM rate-","shift in DR- and BM rate+",
                                         "shift in DR+ and BM rate-","shift in DR+ and BM rate+"),pch=16,
                      col=colo[c(1,2,4,3)],bty="n",pt.cex=1.7,x.intersp=0.7,cex=.8,plot=FALSE)
          -0.8*leg$rect$h/(max(plotinfo$yy)-min(plotinfo$yy))->ii
          legend("topleft",legend=c("shift in DR- and BM rate-","shift in DR- and BM rate+",
                                    "shift in DR+ and BM rate-","shift in DR+ and BM rate+"),pch=16,
                 col=colo[c(1,2,4,3)],bty="n",pt.cex=1.7,x.intersp=0.7,cex=.8,inset=c(0,ii),xpd=TRUE)
        }else{
          if(Ntip(tree)>100) plot(tree, show.tip.label = FALSE) else plot(tree, cex=.8)
        }
        dev.off()
      }
    } else {
      ##### Node #####
      DR->DRrates
      makeL1(tree)->L1
      nod<-list()
      for(k in 1:length(node)){
        as.numeric(names(L1[,which(colnames(L1)==node[k])][which(L1[,which(colnames(L1)==node[k])]!=0)]))[-1]->des
        ldes<-array()
        for(i in 1:length(des)) length(which(des%in%getMommy(tree,des[i])))->ldes[i]
        des[which(ldes<=1)]->des1
        c(getMommy(tree,node[k])[1:2],node[k],des1)->nod[[k]]
      }
      names(nod)<-node
      stack(nod)->paths
      
      unlist(nod)->nods
      if(any(duplicated(nods))) nods[-which(duplicated(nods))]->nods
      if((Ntip(tree)+1)%in%nods) nods[-which(nods==(Ntip(tree)+1))]->nods
      if(any(is.na(nods))) nods[-which(is.na(nods))]->nods
      
      #### DR rates ####
      leaf2N.diff <- array()
      p.single <- array()
      for (j in 1:length(nods)) {
        Ctips <- tips(tree, as.numeric(nods[j]))
        leaf.rates <- DRrates[match(Ctips, rownames(DRrates)),]
        NCrates <- DRrates[-match(names(leaf.rates), rownames(DRrates)),]
        leafR <- mean(abs(leaf.rates))
        NCR <- mean(abs(NCrates))
        leaf2N.diff[j] <- leafR - NCR
        NC <- length(DRrates) - length(leaf.rates)
        C <- length(leaf.rates)
        ran.diffM <- array()
        for (i in 1:nrep) {
          ran.diffM[i] <- mean(sample(abs(DRrates), C))-mean(sample(abs(DRrates),NC))
        }
        p.single[j] <- rank(c(leaf2N.diff[j], ran.diffM[1:(nrep-1)]))[1]/nrep
      }
      names(leaf2N.diff) <- nods
      names(p.single) <- nods
      
      p.single->p.all
      if(!is.null(RR)){ 
        #### RR rates ####
        p.single->pDR.single
        leaf2N.diff->leaf2N.DR.diff
        
        leaf2N.diff <- array()
        p.single <- array()
        for (j in 1:length(nods)) {
          Cleaf <- getDescendants(tree, nods[j])
          Cleaf[which(Cleaf<=Ntip(tree))]<-tree$tip.label[Cleaf[which(Cleaf<=Ntip(tree))]]
          leaf.rates <- rates[match(Cleaf, rownames(rates)),]
          NCrates <- rates[-match(names(leaf.rates), rownames(rates))]
          leaf2N.diff[j] <- mean(abs(leaf.rates)) - mean(abs(NCrates))
          NC <- length(NCrates)
          C <- length(leaf.rates)
          ran.diffM <- array()
          for (i in 1:nrep) {
            ran.diffM[i] <- mean(sample(abs(rates), C)) - mean(sample(abs(rates),
                                                                      NC))
          }
          p.single[j] <- rank(c(leaf2N.diff[j], ran.diffM[1:(nrep -
                                                               1)]))[1]/nrep
        }
        names(leaf2N.diff) <- nods
        names(p.single) <- nods
        
        data.frame("DR.diff"=leaf2N.DR.diff,"p.DR"=pDR.single,"phen.diff"=leaf2N.diff,"p.phen"=p.single)->allres
        
        #if(length(y)>Ntip(tree)) fastBM(tree,nsim=ncol(y))->yB else fastBM(tree)->yB
        sigrat<-array()
        #sigratB<-array()
        for(k in 1:nrow(allres)){
          if(length(y)>Ntip(tree)) {
            cat<-rep("a",nrow(y))
            names(cat)<-rownames(y) 
          }else{
            cat<-rep("a",length(y))
            names(cat)<-names(y)
          }
          cat[tips(tree,as.numeric(rownames(allres)[k]))]<-"nod"
          compare.evol.rates(cbind(y,y),tree,cat,print.progress = FALSE)->com
          #compare.evol.rates(as.matrix(yB),tree,cat,print.progress = FALSE)->comB
          com[[1]]->sigrat[k]
          #comB[[1]]->sigratB[k]
        }
        data.frame(allres,"sigma.ratio"=sigrat)->allres
        #max(sigratB)*1.5->Bsig
        
        #### selected nodes ####
        allres[match(node,rownames(allres)),]->selres
        selres[,2]->psel
        selres[,1]->leaf2N.diff.sel
        names(psel)<-names(leaf2N.diff.sel)<-rownames(selres)
        
        if(length(y)>Ntip(tree)) fastBM(tree,nsim=ncol(y))->yB else fastBM(tree)->yB
        if(length(psel)>1){
          
          sigratB<-array()
          psel1<-array()
          for (j in 1:length(psel)) {
            
            ##### NEW BR NODE
            if(length(y)>Ntip(tree)) {
              cat<-rep("a",nrow(y))
              names(cat)<-rownames(y) 
            }else{
              cat<-rep("a",length(y))
              names(cat)<-names(y)
            }
            cat[tips(tree,as.numeric(names(psel)[j]))]<-"nod"
            compare.evol.rates(cbind(yB,yB),tree,cat,print.progress = FALSE)->comB
            comB[[1]]->sigratB[j]
            ####
            
            unlist(lapply(as.numeric(names(psel)[-j]),function(x) tips(tree,x)))->rem.sel1
            Ctips.sel1<- tips(tree, as.numeric(names(psel)[j]))
            leaf.DR.sel <- DRrates[match(Ctips.sel1, rownames(DRrates)),]
            NCDR.sel <- DRrates[-c(match(names(leaf.DR.sel), rownames(DRrates)),
                                   match(rem.sel1, rownames(DRrates))),]
            leaf2N.DR.sel1 <- mean(abs(leaf.DR.sel)) - mean(abs(NCDR.sel))
            NCDR.sel <- length(NCDR.sel)
            CDR.sel <- length(leaf.DR.sel)
            
            unlist(lapply(as.numeric(names(psel)[-j]),function(x) getDescendants(tree,x)))->rem.sel2
            rem.sel2[which(rem.sel2<=Ntip(tree))]<-tree$tip.label[rem.sel2[which(rem.sel2<=Ntip(tree))]]
            Ctips.sel2 <- getDescendants(tree, as.numeric(names(dwup)[j]))
            Ctips.sel2[which(Ctips.sel2<=Ntip(tree))]<-tree$tip.label[Ctips.sel2[which(Ctips.sel2<=Ntip(tree))]]
            leaf.rates.sel <- rates[match(Ctips.sel2, rownames(rates)),]
            NCrates.sel <- rates[-c(match(names(leaf.rates.sel), rownames(rates)),
                                    match(rem.sel2, rownames(rates))),]
            
            leaf2N.diff.sel <- mean(abs(leaf.rates.sel))-mean(abs(NCrates.sel))
            NC.sel <- length(NCrates.sel)
            C.sel <- length(leaf.rates.sel)
            
            
            ran.DR.sel <- array()
            for (i in 1:nrep) {
              ran.DR.sel[i] <- mean(sample(abs(DRrates), CDR.sel)) - mean(sample(abs(DRrates),NCDR.sel))
              ran.diffM.sel[i] <- mean(sample(abs(rates), C.sel)) - mean(sample(abs(rates),NC.sel))
            }
            psel1[j] <- rank(c(leaf2N.DR.sel1, ran.DR.sel[1:(nrep -1)]))[1]/nrep
            psel2[j] <- rank(c(leaf2N.diff.sel, ran.diffM.sel[1:(nrep -1)]))[1]/nrep
          }
          names(psel1)<-names(psel2)<-names(psel)
          
          data.frame(selres,
                     p.DR.single=psel1[match(rownames(selres),names(psel1))],
                     p.phen.single=psel2[match(rownames(selres),names(psel2))])->selres
        }else{
          
          if(length(y)>Ntip(tree)) {
            cat<-rep("a",nrow(y))
            names(cat)<-rownames(y) 
          }else{
            cat<-rep("a",length(y))
            names(cat)<-names(y)
          }
          cat[tips(tree,as.numeric(names(psel)))]<-"nod"
          compare.evol.rates(cbind(yB,yB),tree,cat,print.progress = FALSE)->comB
          comB[[1]]->sigratB
        }
        names(sigratB)<-names(psel)
        sigratB*2->Bsig
        data.frame(paths,Bsig=Bsig[match(paths[,2],names(Bsig))])->paths

        #### path nodes ####
        allres[-match(node,rownames(allres)),]->allres
        
        allres[which(allres[,2]<=0.025|allres[,2]>=0.975),]->rescut
        if(nrow(rescut)>0) rescut[which(rescut[,4]<=0.025|rescut[,4]>=0.975),]->signres else rescut->signres
        #if(nrow(signres)>0) signres[which(signres[,5]>=Bsig),]->single else signres->single
        data.frame(signres,Bsig=paths[match(rownames(signres),paths[,1]),3])->signres
        if(nrow(signres)>0) signres[which(signres[,5]>=signres[,6]),-6]->single else signres->single
        
        
        if(nrow(single)>1){
          single[order(single[,3]),]->single
          ups.all<-ups <- single[which(single[,2] >= 0.975),2]
          names(ups)<-rownames(single[which(single[,2] >= 0.975),])
          dws.all<-dws <- single[which(single[,2] <= 0.025),2]
          names(dws)<-rownames(single[which(single[,2] <= 0.025),])
          
          allres.diff <- single[,3]
          names(allres.diff)<-rownames(single)
          
          if (length(dws)!=0) {
            s = 1
            repeat{
              d <- which(names(dws) %in% getDescendants(tree, names(dws)[s]))
              if (length(d) > 0) {
                allres.diff[c(match(names(dws[d]),names(allres.diff)),match(names(dws[s]),names(allres.diff)))]->cla
                names(which.max(abs(cla)))->IN
                dws[-match(names(cla[which(names(cla)!=IN)]),names(dws))]->dws
                s=1
              } else {
                dws <- dws
                s = s+1
              }
              if (s > length(dws))  break
            }
          }else dws<-NULL
          
          if (length(ups)!=0) {
            z = 1
            repeat{
              d <- which(names(ups) %in% getDescendants(tree, names(ups)[z]))
              if (length(d) > 0) {
                allres.diff[c(match(names(ups[d]),names(allres.diff)),match(names(ups[z]),names(allres.diff)))]->cla
                names(which.max(abs(cla)))->IN
                ups[-match(names(cla[which(names(cla)!=IN)]),names(ups))]->ups
                z=1
              }else{
                ups <- ups
                z = z+1
              }
              if (z > length(ups))  break
            }
          } else ups<-NULL 
          
          if(length(ups)>=1&length(dws)>=1){
            
            p.single.dws<-array()
            for (j in 1:length(dws)) {
              unlist(lapply(as.numeric(names(ups)),function(x) tips(tree,x)))->rem.ups
              Ctips.dws <- tips(tree, as.numeric(names(dws)[j]))
              leaf.rates.dws <- DRrates[match(Ctips.dws, rownames(DRrates)),]
              
              c(getMommy(tree,names(dws)[j])[1],getDescendants(tree,names(dws)[j])[1])->imm
              if(any(imm%in%names(dws.all))) imm[which(imm%in%names(dws.all))]->bar else bar<-NULL
              if(length(bar)>0){
                roll<-list(); for(w in 1:length(bar)) tips(tree,as.numeric(bar[w]))->roll[[w]]
                unique(c(unlist(roll),names(leaf.rates.dws)))->tm
                NCrates.dws <- DRrates[-c(match(tm, rownames(DRrates)),
                                          match(rem.ups, rownames(DRrates))),]
              }else{
                NCrates.dws <- DRrates[-c(match(names(leaf.rates.dws), rownames(DRrates)),
                                          match(rem.ups, rownames(DRrates))),]
              }
              
              leaf2N.diff.dws <- mean(abs(leaf.rates.dws)) - mean(abs(NCrates.dws))
              NC.dws <- length(NCrates.dws)
              C.dws <- length(leaf.rates.dws)
              
              ran.diffM.dws <- array()
              for (i in 1:nrep) {
                ran.diffM.dws[i] <- mean(sample(abs(DRrates), C.dws)) - mean(sample(abs(DRrates),NC.dws))
              }
              p.single.dws[j] <- rank(c(leaf2N.diff.dws, ran.diffM.dws[1:(nrep -1)]))[1]/nrep
            }
            
            p.single.ups <- array()
            for (j in 1:length(ups)) {
              unlist(lapply(as.numeric(names(dws)),function(x) tips(tree,x)))->rem.dws
              Ctips.ups <- tips(tree, as.numeric(names(ups)[j]))
              leaf.rates.ups <- DRrates[match(Ctips.ups, rownames(DRrates)),]
              
              c(getMommy(tree,names(ups)[j])[1],getDescendants(tree,names(ups)[j])[1])->imm
              if(any(imm%in%names(ups.all))) imm[which(imm%in%names(ups.all))]->bar else bar<-NULL
              if(length(bar)>0){
                roll<-list(); for(w in 1:length(bar)) tips(tree,as.numeric(bar[w]))->roll[[w]]
                unique(c(unlist(roll),names(leaf.rates.ups)))->tm
                NCrates.ups <- DRrates[-c(match(tm, rownames(DRrates)),
                                          match(rem.dws, rownames(DRrates))),]
              }else{
                NCrates.ups <- DRrates[-c(match(names(leaf.rates.ups), rownames(DRrates)),
                                          match(rem.dws, rownames(DRrates))),]
              }
              
              
              leaf2N.diff.ups <- mean(abs(leaf.rates.ups)) - mean(abs(NCrates.ups))
              NC.ups <- length(NCrates.ups)
              C.ups <- length(leaf.rates.ups)
              
              ran.diffM.ups <- array()
              for (i in 1:nrep) {
                ran.diffM.ups[i] <- mean(sample(abs(DRrates), C.ups)) - mean(sample(abs(DRrates),NC.ups))
              }
              p.single.ups[j] <- rank(c(leaf2N.diff.ups, ran.diffM.ups[1:(nrep -1)]))[1]/nrep
            }
            
            if(!all(c(p.single.dws>0.05,p.single.ups<0.95))){
              if(any(p.single.dws>0.05)) dws[-which(p.single.dws>0.05)]->dws
              if(any(p.single.ups<0.95)) ups[-which(p.single.ups<0.95)]->ups
            }
          }
          
          c(dws,ups)->dwup
          p.single.dwup<-array()
          for (j in 1:length(dwup)) {
            unlist(lapply(as.numeric(names(dwup)[-j]),function(x) getDescendants(tree,x)))->rem
            rem[which(rem<=Ntip(tree))]<-tree$tip.label[rem[which(rem<=Ntip(tree))]]
            Ctips.foc <- getDescendants(tree, as.numeric(names(dwup)[j]))
            Ctips.foc[which(Ctips.foc<=Ntip(tree))]<-tree$tip.label[Ctips.foc[which(Ctips.foc<=Ntip(tree))]]
            leaf.rates.foc <- rates[match(Ctips.foc, rownames(rates)),]
            NCrates.foc <- rates[-c(match(names(leaf.rates.foc), rownames(rates)),
                                    match(rem, rownames(rates))),]
            
            leaf2N.diff.foc <- mean(abs(leaf.rates.foc))-mean(abs(NCrates.foc))
            NC.foc <- length(NCrates.foc)
            C.foc <- length(leaf.rates.foc)
            
            ran.diffM.foc <- array()
            for (i in 1:nrep) {
              ran.diffM.foc[i] <- mean(sample(abs(rates), C.foc)) - mean(sample(abs(rates),NC.foc))
            }
            p.single.dwup[j] <- rank(c(leaf2N.diff.foc, ran.diffM.foc[1:(nrep -1)]))[1]/nrep
          }
          names(p.single.dwup)<-names(dwup)
          
          if(any(p.single.dwup>=0.975|p.single.dwup<=0.025)) 
            single<-single[which(rownames(single)%in%names(which(p.single.dwup>=0.975|p.single.dwup<=0.025))),] else single<-NULL
          
        }
        
        if(!is.null(single)) if(nrow(single)<1) single<-NULL
        
        list(selres,single,Bsig)->shift.res
        names(shift.res)<-c("focal.nodes","path.nodes","max.s2ratio.bm")
        
      }else{ #### BM rates ####
        #if(length(y)>Ntip(tree)) fastBM(tree,nsim=ncol(y))->yB else fastBM(tree)->yB
        aovres<-list()
        #sigratB<-array()
        for(k in 1:length(p.all)){
          if(length(y)>Ntip(tree)) {
            cat<-rep("a",nrow(y))
            names(cat)<-rownames(y) 
          }else{
            cat<-rep("a",length(y))
            names(cat)<-names(y)
          }
          cat[tips(tree,as.numeric(names(p.all[k])))]<-"nod"
          compare.evol.rates(cbind(y,y),tree,cat,print.progress = FALSE)->com
          c(com[[5]][which(names(com[[5]])=="nod")],com[[5]][-which(names(com[[5]])=="nod")])->sigma.g
          c(unlist(com[1:2]),sigma.g)->aovres[[k]]
          #compare.evol.rates(as.matrix(yB),tree,cat,print.progress = FALSE)[[1]]->sigratB[k]
        }
        do.call(rbind,aovres)->aovres
        rownames(aovres)<-names(p.all)
        colnames(aovres)<-c("sigma.ratio","p.sigma","sigma.node","sigma.back")
        #max(sigratB)*1.5->Bsig
        
        data.frame("DR.diff"=leaf2N.diff,"p.DR"=p.all,aovres)->allres
        
        #### selected nodes ####
        allres[match(node,rownames(allres)),]->selres
        p.single[match(node,names(p.single))]->psel
        leaf2N.diff[match(node,names(leaf2N.diff))]->leaf2N.diff.sel
        if(length(y)>Ntip(tree)) fastBM(tree,nsim=ncol(y))->yB else fastBM(tree)->yB
        
        if(length(psel)>1){
          psel1<-array()
          sigratB<-array()
          for (j in 1:length(psel)) {
            
            ##### NEW BR NODE
            if(length(y)>Ntip(tree)) {
              cat<-rep("a",nrow(y))
              names(cat)<-rownames(y) 
            }else{
              cat<-rep("a",length(y))
              names(cat)<-names(y)
            }
            cat[tips(tree,as.numeric(names(psel)[j]))]<-"nod"
            compare.evol.rates(cbind(yB,yB),tree,cat,print.progress = FALSE)->comB
            comB[[1]]->sigratB[j]
            
            unlist(lapply(as.numeric(names(psel)[-j]),function(x) tips(tree,x)))->rem.sel
            Ctips.sel<- tips(tree, as.numeric(names(psel)[j]))
            leaf.rates.sel <- DRrates[match(Ctips.sel, rownames(DRrates)),]
            NCrates.sel <- DRrates[-c(match(names(leaf.rates.sel), rownames(DRrates)),
                                      match(rem.sel, rownames(DRrates))),]
            leafR.sel <- mean(abs(leaf.rates.sel))
            NCR.sel <- mean(abs(NCrates.sel))
            leaf2N.diff.sel1 <- leafR.sel - NCR.sel
            NC.sel <- length(NCrates.sel)
            C.sel <- length(leaf.rates.sel)
            
            ran.diffM.sel <- array()
            for (i in 1:nrep) {
              ran.diffM.sel[i] <- mean(sample(abs(DRrates), C.sel)) - mean(sample(abs(DRrates),NC.sel))
            }
            psel1[j] <- rank(c(leaf2N.diff.sel1, ran.diffM.sel[1:(nrep -1)]))[1]/nrep
            
          }
          names(psel1)<-names(psel)
          
          lapply(as.numeric(names(psel)),function(x) tips(tree,x))->sta
          names(sta)<-names(psel)
          stack(sta)->sta1
          
          if(length(y)>Ntip(tree)) {
            cat<-rep("back",nrow(y))
            names(cat)<-rownames(y) 
          }else{
            cat<-rep("back",length(y))
            names(cat)<-names(y)
          }
          cat[match(sta1[,1],names(cat))]<-paste("g",sta1[,2],sep="")
          
          compare.evol.rates(cbind(y,y),tree,cat,print.progress = FALSE)->cer
          as.matrix(cer$pairwise.pvalue)[-1,1]->pcer
          gsub("g","",names(pcer))->names(pcer)
          
          data.frame(selres[,c(1:3,5:6,4)],
                     p.DR.single=psel1[match(rownames(selres),names(psel1))],
                     p.sigma.single=pcer[match(rownames(selres),names(pcer))])->selres
        }else{
          
          selres[,c(1:3,5:6,4)]->selres
          
          if(length(y)>Ntip(tree)) {
            cat<-rep("a",nrow(y))
            names(cat)<-rownames(y) 
          }else{
            cat<-rep("a",length(y))
            names(cat)<-names(y)
          }
          cat[tips(tree,as.numeric(names(psel)))]<-"nod"
          compare.evol.rates(cbind(yB,yB),tree,cat,print.progress = FALSE)->comB
          comB[[1]]->sigratB
        }
        names(sigratB)<-names(psel)
        sigratB*2->Bsig  
        
        #### path nodes ####
        allres[-match(node,rownames(allres)),]->allres
        p.single[-match(node,names(p.single))]->p.single
        
        data.frame(paths,Bsig=Bsig[match(paths[,2],names(Bsig))])->paths
        data.frame(allres,Bsig=paths[match(rownames(allres),paths[,1]),3])->allres
        
        if(length(which(allres[,4]<=0.05))==0)  p.single<-NULL else p.single[which(allres[,4]<=0.05&allres[,3]>allres[,7])]->p.single
        
        if (length(p.single[p.single >= 0.975 | p.single <=0.025])==0)  p.single<-NULL
        
        if (length(p.single[p.single >= 0.975 | p.single <=
                            0.025]) == 1) {
          p.single <- p.single[p.single >= 0.975 | p.single <=
                                 0.025]
          leaf2N.diff <- leaf2N.diff[match(names(p.single),
                                           names(leaf2N.diff))]
        }
        
        if (length(p.single[p.single >= 0.975 | p.single <=
                            0.025]) >= 2)  {
          
          p.single <- p.single[p.single >= 0.975 | p.single <= 0.025]
          allres.diff <- allres[match(names(p.single), rownames(allres)),3]
          names(allres.diff)<-names(p.single)
          
          ups <- p.single[p.single >= 0.975]
          dws <- p.single[p.single <= 0.025]
          
          ups.all <- ups <- ups[na.omit(match(names(allres.diff[order(allres.diff,
                                                                      decreasing = FALSE)]), names(ups)))]
          dws.all <- dws <- dws[na.omit(match(names(allres.diff[order(allres.diff,
                                                                      decreasing = FALSE)]), names(dws)))]
          
          if (length(dws)!=0) {
            s = 1
            repeat{
              d <- which(names(dws) %in% getDescendants(tree, names(dws)[s]))
              if (length(d) > 0) {
                allres.diff[c(match(names(dws[d]),names(allres.diff)),match(names(dws[s]),names(allres.diff)))]->cla
                names(which.max(abs(cla)))->IN
                dws[-match(names(cla[which(names(cla)!=IN)]),names(dws))]->dws
                s=1
              }else{
                dws <- dws
                s = s+1
              }
              if (s > length(dws))  break
            }
          } else dws<-NULL
          if (length(ups)!=0) {
            z = 1
            repeat{
              d <- which(names(ups) %in% getDescendants(tree, names(ups)[z]))
              if (length(d) > 0) {
                allres.diff[c(match(names(ups[d]),names(allres.diff)),match(names(ups[z]),names(allres.diff)))]->cla
                names(which.max(abs(allres.diff[c(match(names(ups[d]),names(allres.diff)),match(names(ups[z]),names(allres.diff)))])))->IN
                ups[-match(names(cla[which(names(cla)!=IN)]),names(ups))]->ups
                z=1
              }else{
                ups <- ups
                z = z+1
              }
              if (z > length(ups))  break
            }
          } else ups<-NULL
          
          if(length(ups)>=1&length(dws)>=1){
            p.single.dws<-array()
            for (j in 1:length(dws)) {
              unlist(lapply(as.numeric(names(ups)),function(x) tips(tree,x)))->rem.ups
              Ctips.dws <- tips(tree, as.numeric(names(dws)[j]))
              leaf.rates.dws <- DRrates[match(Ctips.dws, rownames(DRrates)),]
              
              c(getMommy(tree,names(dws)[j])[1],getDescendants(tree,names(dws)[j])[1])->imm
              if(any(imm%in%names(dws.all))) imm[which(imm%in%names(dws.all))]->bar else bar<-NULL
              if(length(bar)>0){
                roll<-list(); for(w in 1:length(bar)) tips(tree,as.numeric(bar[w]))->roll[[w]]
                unique(c(unlist(roll),names(leaf.rates.dws)))->tm
                NCrates.dws <- DRrates[-c(match(tm, rownames(DRrates)),
                                          match(rem.ups, rownames(DRrates))),]
              }else{
                NCrates.dws <- DRrates[-c(match(names(leaf.rates.dws), rownames(DRrates)),
                                          match(rem.ups, rownames(DRrates))),]
              }
              
              leafR.dws <- mean(abs(leaf.rates.dws))
              NCR.dws <- mean(abs(NCrates.dws))
              leaf2N.diff.dws <- leafR.dws - NCR.dws
              NC.dws <- length(NCrates.dws)
              C.dws <- length(leaf.rates.dws)
              
              ran.diffM.dws <- array()
              for (i in 1:nrep) {
                ran.diffM.dws[i] <- mean(sample(abs(DRrates), C.dws)) - mean(sample(abs(DRrates),NC.dws))
              }
              p.single.dws[j] <- rank(c(leaf2N.diff.dws[j], ran.diffM.dws[1:(nrep -1)]))[1]/nrep
            }
            
            p.single.ups <- array()
            for (j in 1:length(ups)) {
              unlist(lapply(as.numeric(names(dws)),function(x) tips(tree,x)))->rem.dws
              Ctips.ups <- tips(tree, as.numeric(names(ups)[j]))
              leaf.rates.ups <- DRrates[match(Ctips.ups, rownames(DRrates)),]
              
              c(getMommy(tree,names(ups)[j])[1],getDescendants(tree,names(ups)[j])[1])->imm
              if(any(imm%in%names(ups.all))) imm[which(imm%in%names(ups.all))]->bar  else bar<-NULL
              if(length(bar)>0){
                roll<-list(); for(w in 1:length(bar)) tips(tree,as.numeric(bar[w]))->roll[[w]]
                unique(c(unlist(roll),names(leaf.rates.ups)))->tm
                NCrates.ups <- DRrates[-c(match(tm, rownames(DRrates)),
                                          match(rem.dws, rownames(DRrates))),]
              }else{
                NCrates.ups <- DRrates[-c(match(names(leaf.rates.ups), rownames(DRrates)),
                                          match(rem.dws, rownames(DRrates))),]
              }
              
              leafR.ups <- mean(abs(leaf.rates.ups))
              NCR.ups <- mean(abs(NCrates.ups))
              leaf2N.diff.ups <- leafR.ups - NCR.ups
              NC.ups <- length(NCrates.ups)
              C.ups <- length(leaf.rates.ups)
              
              ran.diffM.ups <- array()
              for (i in 1:nrep) {
                ran.diffM.ups[i] <- mean(sample(abs(DRrates), C.ups)) - mean(sample(abs(DRrates),NC.ups))
              }
              p.single.ups[j] <- rank(c(leaf2N.diff.ups[j], ran.diffM.ups[1:(nrep -1)]))[1]/nrep
            }
            
            if(!all(c(p.single.dws>0.05,p.single.ups<0.95))){
              if(any(p.single.dws>0.05)) dws[-which(p.single.dws>0.05)]->dws
              if(any(p.single.ups<0.95)) ups[-which(p.single.ups<0.95)]->ups
            }
            
          }
          
          p.single <- p.single[which(names(p.single) %in%
                                       names(c(ups, dws)))]
          
          lapply(names(p.single),function(x) tips(tree,as.numeric(x)))->sta
          names(sta)<-names(p.single)
          stack(sta)->sta1
          
          if(length(y)>Ntip(tree)) {
            cat<-rep("back",nrow(y))
            names(cat)<-rownames(y) 
          }else{
            cat<-rep("back",length(y))
            names(cat)<-names(y)
          }
          cat[match(sta1[,1],names(cat))]<-paste("g",sta1[,2],sep="")
          
          compare.evol.rates(cbind(y,y),tree,cat,print.progress = FALSE)->cer
          if(inherits(cer$pairwise.pvalue,"dist")) as.matrix(cer$pairwise.pvalue)[-1,1]->pcer else {
            cer$pairwise.pvalue->pcer
            names(pcer)<-cer$groups[2]
          }
          if(length(which(pcer<=0.05))>0)
            p.single <- p.single[match(gsub("g","",names(which(pcer<=0.05))),names(p.single))] else p.single<-NULL
          
          leaf2N.diff <- leaf2N.diff[match(names(p.single),
                                           names(leaf2N.diff))]
          
        }
        
        if(length(p.single)>1){
          data.frame("DR.diff"=leaf2N.diff[match(names(p.single),names(leaf2N.diff))],
                     "p.DR"=p.single,aovres[match(names(p.single),rownames(aovres)),c(1,3,4,2)])->single.clades
        }else{
          if(is.null(p.single)) single.clades<-NULL else {
            as.data.frame(matrix(c(unname(leaf2N.diff[match(names(p.single),names(leaf2N.diff))]),
                                   unname(p.single),aovres[match(names(p.single),rownames(aovres)),c(1,3,4,2)]),nrow=1))->single.clades
            
            rownames(single.clades)<-names(p.single)
            colnames(single.clades)<-c("DR.diff","p.DR",colnames(aovres)[c(1,3,4,2)])
          }
        }
        
        list(selres,single.clades,Bsig)->shift.res
        names(shift.res)<-c("focal.nodes","path.nodes","max.s2ratio.bm")
      }
    }
  } else shift.res<-NULL
  
  #### Trend ####
  if(is.null(DR.log)){
    if(shapiro.test(DR)$p<0.05){ 
      log(DR)->DR 
      DR.log<-TRUE
    }else DR.log<-FALSE 
  }else{
    if(isTRUE(DR.log)) log(DR)->DR
  }
  
  if(select.test=="both" | select.test=="trend"){
    pdf(file=paste(foldername, "Character associated diversification plot.pdf",sep="/"))
    if(length(y)>Ntip(tree)) {
      if(ncol(y)<7) par(mfrow=c(ceiling(ncol(y)/2),2),mar = c(3.5, 3.5, 1.5, 1)) else par(mfrow=c(3,2),mar = c(3.5, 3.5, 1.5, 1))
      if(!is.null(colnames(y))) ynams<-colnames(y) else ynams<-paste("y",seq(1:ncol(y)),sep="")
    }
    trend.res<-list()
    for(w in 1:ncol(as.matrix(y))){
      RRphylo:::range01(as.matrix(y)[,w])->yR
      fastBM(tree,nsim=nrep)->yBB
      coef(summary(lm(yR~DR)))[2]->corr
      predict(lm(yR~DR))[which.min(DR)]->miny
      apply(yBB,2,function(x) RRphylo:::range01(x))->yBB
      apply(yBB,2,function(x) coef(summary(lm(x~DR)))[2,1])->corBB
      rank(c(corr[1],corBB)[1:nrep])[1]/nrep->p
      as.data.frame(cbind(slope=coef(summary(lm(yR~DR)))[2],
                          p.real=coef(summary(lm(yR~DR)))[2,4],
                          p.random=p))->trend.res[[w]]
      
      sapply(quantile(corBB,c(0.025,0.975)), function(x) corBB[which.min(abs(corBB-x))])->quaBB
      data.frame(min=(miny-quaBB[1]*min(DR))+quaBB[1]*DR,max=(miny-quaBB[2]*min(DR))+quaBB[2]*DR,x=DR)->CI
      CI[order(CI$x),]->CI
      
      if(isTRUE(DR.log)) xxlab="log DR" else xxlab="DR"
      plot(DR,yR,mgp = c(2, 0.5, 0),xlab="",ylab="phenotype",cex.lab=1.2)
      if(length(y)>Ntip(tree)) title(main=ynams[w],xlab=xxlab,mgp=c(1.5,0.5,0),cex.lab=1.2) else title(xlab=xxlab,mgp=c(1.5,0.5,0),cex.lab=1.2)
      polygon(c(CI[,3],rev(CI[,3])),c(CI[,1],rev(CI[,2])),border=NA,col=rgb(0.5,0.5,0.5,0.4))
      ablineclip(lm(yR~DR),x1=min(DR),x2=max(DR),col="aquamarine",lwd=3,lend=0)
      
    }
    dev.off()
    if(length(y)>Ntip(tree)){
      names(trend.res)<-ynams
      do.call(rbind,trend.res)->trend.res
    }else unlist(trend.res)->trend.res
  }else trend.res<-NULL
  
  list(tree=tree,DR.log=DR.log,shift.res=shift.res,trend.res=trend.res)->res
  return(res)
}
