R/search.trend.R

In RRphylo: Phylogenetic Ridge Regression Methods for Comparative Studies

#'@title Searching for evolutionary trends in phenotypes and rates
#'@description This function searches for evolutionary trends in the phenotypic
#'  mean and the evolutionary rates for the entire tree and individual clades.
#'@usage search.trend(RR,y,x1=NULL,x1.residuals = FALSE,
#'  node=NULL,cov=NULL,nsim=100,clus=0.5,ConfInt=NULL,filename=NULL)
#'@param RR an object produced by \code{\link{RRphylo}}.
#'@param y the named vector (or matrix if multivariate) of phenotypes.
#'@param x1 the additional predictor to be specified if the RR object has been
#'  created using an additional predictor (i.e. multiple version of
#'  \code{RRphylo}). \code{'x1'} vector must be as long as the number of nodes
#'  plus the number of tips of the tree, which can be obtained by running
#'  \code{RRphylo} on the predictor as well, and taking the vector of ancestral
#'  states and tip values to form the \code{x1}. Note: only one predictor at
#'  once can be specified.
#'@param x1.residuals logical specifying whether the residuals of regression
#'  between \code{y} and \code{x1} should be inspected for a phenotypic trend
#'  (see details and examples below). Default is \code{FALSE}.
#'@param node the node number of individual clades to be specifically tested and
#'  contrasted to each other. It is \code{NULL} by default. Notice the node
#'  number must refer to the dichotomic version of the original tree, as
#'  produced by \code{RRphylo}.
#'@param cov the covariate values to be specified if the RR object has been
#'  created using a  covariate for rates calculation.  As for \code{RRphylo},
#'  \code{'cov'} must be as long as the number of nodes plus the number of tips
#'  of the tree, which can be obtained by running \code{RRphylo} on the
#'  covariate as well, and taking the vector of ancestral states and tip values
#'  to form the covariate (see the example below).
#'@param nsim number of simulations to be performed. It is set at 100 by
#'  default.
#'@param clus the proportion of clusters to be used in parallel computing. To
#'  run the single-threaded version of \code{search.trend} set \code{clus} = 0.
#'@param ConfInt is deprecated.
#'@param filename is deprecated. \code{search.trend} does not return plots
#'  anymore, check the function \code{\link{plotTrend}} instead.
#'@return The function returns a list object containing:
#'@return \strong{$trends.data} a 'RRphyloList' object including:
#'  \enumerate{\item{\code{$phenotypeVStime}}: a data frame of phenotypic values
#'  (or \code{y} versus \code{x1} regression residuals if
#'  \code{x1.residuals=TRUE}) and their distance from the tree root for each
#'  node (i.e. ancestral states) and tip of the tree.
#'  \item{\code{$absrateVStime}}: a data frame of \code{RRphylo} rates and the
#'  distance from the tree root (age). If y is multivariate, it also includes
#'  the multiple rates for each y vector. If \code{node} is specified, each
#'  branch is classified as belonging to an indicated clade.
#'  \item{\code{$rescaledrateVStime}}: a data frame of rescaled \code{RRphylo}
#'  rates and the distance from the tree root (age). If y is multivariate, it
#'  also includes the multiple rates for each y vector. If \code{node} is
#'  specified, each branch is classified as belonging to an indicated clade. NAs
#'  correspond either to very small values or to outliers which are excluded
#'  from the analysis.}
#'@return \strong{$phenotypic.regression} results of phenotype (\code{y} versus
#'  \code{x1} regression residuals) versus age regression. It reports a p-value
#'  for the regression slope between the variables (p.real), a p-value computed
#'  contrasting the real slope to Brownian motion simulations (p.random), and a
#'  parameter indicating the deviation of the phenotypic mean from the root
#'  value in terms of the number of standard deviations of the trait
#'  distribution (dev). dev is 0 under Brownian Motion. Only p.random should be
#'  inspected to assess significance.
#'@return \strong{$rate.regression} results of the rates (rescaled absolute
#'  values) versus age regression. It reports a p-value for the regression
#'  between the variables (p.real), a p-value computed contrasting the real
#'  slope to Brownian motion simulations (p.random), and a parameter indicating
#'  the ratio between the range of phenotypic values and the range of such
#'  values halfway along the tree height, divided to the same figure under
#'  Brownian motion (spread). spread is 1 under Brownian Motion. Only p.random
#'  should be inspected to assess significance.
#'@return \strong{$ConfInts} a 'RRphyloList' object including the 95\%
#'  confidence intervals around regression slopes of phenotypes and rates (both
#'  rescaled and unscaled absolute rates) produced according to the Brownian
#'  motion model of evolution.
#'@return If specified, individual nodes are tested as the whole tree, the
#'  results are summarized in the objects:
#'@return \strong{$node.phenotypic.regression} results of phenotype (or \code{y}
#'  versus \code{x1} regression residuals) versus age regression through node.
#'  It reports the slope for the regression between the variables at node
#'  (slope), a p-value computed contrasting the real slope to Brownian motion
#'  simulations (p.random), the difference between estimated marginal means
#'  predictions for the group and for the rest of the tree (emm.difference), and
#'  a p-value for the emm.difference (p.emm).
#'@return \strong{$node.rate.regression} results of the rates (absolute values)
#'  versus age regression through node. It reports the difference between
#'  estimated marginal means predictions for the group and for the rest of the
#'  tree (emm.difference), a p-value for the emm.difference (p.emm), the
#'  regression slopes for the group (slope.node) and for the rest of the tree
#'  (slope.others), and a p-value for the difference between such slopes
#'  (p.slope).
#'@return If more than one node is specified, the object
#'  \strong{$group.comparison} reports the same results as
#'  $node.phenotypic.regression and $node.rate.regression obtained by comparing
#'  individual clades to each other.
#'@author Silvia Castiglione, Carmela Serio, Pasquale Raia, Alessandro
#'  Mondanaro, Marina Melchionna, Mirko Di Febbraro, Antonio Profico, Francesco
#'  Carotenuto
#'@details The function simultaneously returns the regression of phenotypes and
#'  phenotypic evolutionary rates against age tested against Brownian motion
#'  simulations to assess significance. To this aim rates are rescaled in the
#'  0-1 range and then logged. To assess significance, slopes are compared to a
#'  family of simulated slopes (BMslopes, where the number of simulations is
#'  equal to \code{nsim}), generated under the Brownian motion, using the
#'  \code{fastBM} function in the package \pkg{phytools}. Individual nodes are
#'  compared to the rest of the tree in different ways depending on whether
#'  phenotypes or rates (always unscaled in this case) versus age regressions
#'  are tested. With the former, the regression slopes for individual clades and
#'  the slope difference between clades is contrasted to slopes obtained through
#'  Brownian motion simulations. For the latter, regression models are tested
#'  and contrasted to each other referring to estimated marginal means, by using
#'  the \code{emmeans} function in the package \pkg{emmeans}.
#'
#'  The \href{../doc/RRphylo.html#predictor}{multiple regression version of
#'  RRphylo} allows to incorporate the effect of an additional predictor in the
#'  computation of evolutionary rates without altering the ancestral character
#'  estimation. Thus, when a multiple \code{RRphylo} output is fed to
#'  \code{search.trend}, the predictor effect is accounted for on the absolute
#'  evolutionary rates, but not on the phenotype. However, in some situations
#'  the user might want to factor out the predictor effect on phenotypes as
#'  well. Under the latter circumstance, by setting the argument
#'  \code{x1.residuals = TRUE}, the residuals of the response to predictor
#'  regression are used as to represent the phenotype.
#'@importFrom stats as.formula coef resid density predict cor
#'@importFrom phytools nodeHeights
#'@importFrom parallel makeCluster detectCores stopCluster
#'@importFrom doParallel registerDoParallel
#'@importFrom foreach foreach %dopar%
#'@importFrom utils combn
#'@importFrom emmeans emmeans emtrends
#'@export
#'@seealso \href{../doc/search.trend.html}{\code{search.trend} vignette}
#'@seealso \href{../doc/Plotting-tools.html}{\code{plotTrend} vignette}
#'@seealso \code{\link{plotTrend}}
#'@references Castiglione, S., Serio, C., Mondanaro, A., Di Febbraro, M.,
#'  Profico, A., Girardi, G., & Raia, P. (2019) Simultaneous detection of
#'  macroevolutionary patterns in phenotypic means and rate of change with and
#'  within phylogenetic trees including extinct species. \emph{PLoS ONE}, 14:
#'  e0210101. https://doi.org/10.1371/journal.pone.0210101
#' @examples
#'  \dontrun{
#' data("DataOrnithodirans")
#' DataOrnithodirans$treedino->treedino
#' DataOrnithodirans$massdino->massdino
#' cc<- 2/parallel::detectCores()
#'
#' # Extract Pterosaurs tree and data
#' library(ape)
#' extract.clade(treedino,746)->treeptero
#' massdino[match(treeptero$tip.label,names(massdino))]->massptero
#' massptero[match(treeptero$tip.label,names(massptero))]->massptero
#'
#' # Case 1. "RRphylo" whitout accounting for the effect of a covariate
#' RRphylo(tree=treeptero,y=log(massptero),clus=cc)->RRptero
#'
#' # Case 1.1. "search.trend" whitout indicating nodes to be tested for trends
#' search.trend(RR=RRptero, y=log(massptero), nsim=100, clus=cc,cov=NULL,node=NULL)
#'
#' # Case 1.2. "search.trend" indicating nodes to be specifically tested for trends
#' search.trend(RR=RRptero, y=log(massptero), nsim=100, node=143, clus=cc,cov=NULL)
#'
#'
#' # Case 2. "RRphylo" accounting for the effect of a covariate
#' # "RRphylo" on the covariate in order to retrieve ancestral state values
#' RRphylo(tree=treeptero,y=log(massptero),clus=cc)->RRptero
#' c(RRptero$aces,log(massptero))->cov.values
#' names(cov.values)<-c(rownames(RRptero$aces),names(massptero))
#' RRphylo(tree=treeptero,y=log(massptero),cov=cov.values,clus=cc)->RRpteroCov
#'
#' # Case 2.1. "search.trend" whitout indicating nodes to be tested for trends
#' search.trend(RR=RRpteroCov, y=log(massptero), nsim=100, clus=cc,cov=cov.values)
#'
#' # Case 2.2. "search.trend" indicating nodes to be specifically tested for trends
#' search.trend(RR=RRpteroCov, y=log(massptero), nsim=100, node=143, clus=cc,cov=cov.values)
#'
#'
#' # Case 3. "search.trend" on multiple "RRphylo"
#' data("DataCetaceans")
#' DataCetaceans$treecet->treecet
#' DataCetaceans$masscet->masscet
#' DataCetaceans$brainmasscet->brainmasscet
#' DataCetaceans$aceMyst->aceMyst
#'
#' drop.tip(treecet,treecet$tip.label[-match(names(brainmasscet),treecet$tip.label)])->treecet.multi
#' masscet[match(treecet.multi$tip.label,names(masscet))]->masscet.multi
#'
#' RRphylo(tree=treecet.multi,y=masscet.multi,clus=cc)->RRmass.multi
#' RRmass.multi$aces[,1]->acemass.multi
#' c(acemass.multi,masscet.multi)->x1.mass
#'
#' RRphylo(tree=treecet.multi,y=brainmasscet,x1=x1.mass,clus=cc)->RRmulti
#'
#' # incorporating the effect of body size at inspecting trends in absolute evolutionary rates
#' search.trend(RR=RRmulti, y=brainmasscet,x1=x1.mass,clus=cc)
#'
#' # incorporating the effect of body size at inspecting trends in both absolute evolutionary
#' # rates and phenotypic values (by using brain versus body mass regression residuals)
#' search.trend(RR=RRmulti, y=brainmasscet,x1=x1.mass,x1.residuals=TRUE,clus=cc)
#'    }


search.trend<-function (RR,y,
                        x1=NULL,x1.residuals=FALSE,
                        node = NULL, cov = NULL,
                        nsim = 100, clus = 0.5, ConfInt = NULL,
                        filename=NULL)
{
  # require(ape)
  # require(phytools)
  # require(stats4)
  # require(foreach)
  # require(doParallel)
  # require(parallel)
  # require(nlme)
  # require(RColorBrewer)
  # require(emmeans)
  # require(outliers)
  # require(car)

  if (!requireNamespace("RColorBrewer", quietly = TRUE)) {
    stop("Package \"RColorBrewer\" needed for this function to work. Please install it.",
         call. = FALSE)
  }

  if (!requireNamespace("car", quietly = TRUE)) {
    stop("Package \"car\" needed for this function to work. Please install it.",
         call. = FALSE)
  }

  if(!missing(filename))
    warning("The argument filename is deprecated. Check the function plotTrend to plot results",immediate. = TRUE)

  if(!missing(ConfInt))
    warning("The argument ConfInt is deprecated.",immediate. = TRUE)


  range01 <- function(x, ...){(x - min(x, ...)) / (max(x, ...) - min(x, ...))}

  t <- RR$tree
  if(min(diag(vcv(t)))/max(diag(vcv(t)))>=0.9) stop("not enough fossil information")
  rates <- RR$rates
  betas <- RR$multiple.rates
  aceRR <- RR$aces
  L <- RR$tip.path
  L1 <- RR$node.path
  if (!is.null(ncol(y))&is.null(rownames(y))) stop("The matrix of phenotypes needs to be named")
  if (is.null(ncol(y))&is.null(names(y))) stop("The vector of phenotypes needs to be named")
  ynam<-deparse(substitute(y))
  # if(is.null(nrow(y))) y <- treedata(t, y, sort = TRUE)[[2]][,1] else y <- treedata(t, y, sort = TRUE)[[2]]
  y <- treedataMatch(t, y)[[1]]

  H <- max(nodeHeights(t))
  eds <- data.frame(leaf = c(Ntip(t)+1,t$edge[, 2]),height=c(0,nodeHeights(t)[, 2]))
  eds[which(eds[,1]<=Ntip(t)),1] <- t$tip.label[eds[which(eds[,1]<=Ntip(t)),1]]

  if (ncol(y) > 1) {
    y.multi <- (L %*% rates)
    y.multi[match(rownames(y),rownames(y.multi)),,drop=FALSE]->y.multi
    aceRR.multi <- (L1 %*% rates[1:Nnode(t), ])
    nodes <- cbind(aceRR,aceRR.multi)
    if(is.null(colnames(y))) colnames(nodes)<- c(paste("y", seq(1,ncol(y)),sep = ""),"y.multi") else
      colnames(nodes)<-c(colnames(y),"y.multi")
    P <- rbind(nodes, cbind(y,y.multi))
    if(isTRUE(x1.residuals)) apply(P,2,function(x) residuals(lm(x~x1)))->P

    rates <- as.data.frame(rates)
    betas <- as.data.frame(betas)

    rate.data <- data.frame(betas = betas[match(eds[, 1], rownames(betas)),],
                            rate = rates[match(eds[, 1], rownames(rates)),],
                            age = eds[, 2])
    if(is.null(colnames(y))) colnames(rate.data)[1:ncol(y)] <- paste("betas", seq(1,ncol(y)), sep = "") else
      colnames(rate.data)[1:ncol(y)]<-colnames(y)

  }else{
    P <- rbind(aceRR, y)
    if(isTRUE(x1.residuals)) as.matrix(residuals(lm(P~x1)))->P
    # colnames(P)<-"y"
    colnames(P)<-ynam

    rate.data <- data.frame(rate = rates[match(eds[, 1], rownames(rates)),],
                            age = eds[, 2])
    rownames(rate.data) <- rownames(rates)[match(eds[, 1], rownames(rates))]

  }

  rate.data[,ncol(rate.data)]+L[1,1]->rate.data[,ncol(rate.data)]
  phen.data <- data.frame(P[match(rownames(rate.data), rownames(P)),,drop=FALSE],age=rate.data$age,check.names = FALSE)
  rate.dataRES<-rate.data
  phen.dataRES <- phen.data
  {##### Rate Trend Real Multi #####
    rate.reg <- scalrat.data<-rate.coef <-list()
    e1<-array()
    for (i in 1:(ncol(rate.data)-1)) {
      bet <- rate.data[, i,drop=FALSE]
      age <- rate.data[, ncol(rate.data),drop=FALSE]
      abs(bet)->rts->rtsA
      log(range01(rts))->rts

      c(which(rts=="-Inf"),which(age=="-Inf"))->outs
      if(length(outs)>0){
        rts[-outs,,drop=FALSE]->rts
        age[-outs,,drop=FALSE]->age
      }
      age->ageC

      sd(range01(rtsA[ageC<0.5*max(ageC)]))/sd(range01(rtsA)[ageC>0.5*max(ageC)])->e1[i]

      if(!is.null(x1)){
        rts[-1,,drop=FALSE]->rts
        age[-1,,drop=FALSE]->age

        car::outlierTest(lm(as.matrix(rts)~as.matrix(age)))->ouT
        if(any(ouT$bonf.p<=0.05)){
          rts[-match(names(which(ouT$bonf.p<=0.05)),rownames(rts)),,drop=FALSE]->rts
          age[-match(names(which(ouT$bonf.p<=0.05)),rownames(age)),,drop=FALSE]->age
        }
      }else{
        lm(as.matrix(rts)~as.matrix(age))->bb

        residuals(bb)[order(residuals(bb),decreasing=TRUE)][1:(Ntip(t)/15)]->resout
        if((Ntip(t)+1)%in%names(resout)){
          rts[-1,,drop=FALSE]->rts
          age[-1,,drop=FALSE]->age
        }
      }
      lm(as.matrix(rts)~as.matrix(age))->regr.1
      rate.coef[[i]] <- coef(summary(regr.1))[2, c(1, 4)]
      rate.reg[[i]] <- regr.1

      scalrat.data[[i]]<-rts[match(rownames(rate.data),rownames(rts)),,drop=FALSE]
      rownames(scalrat.data[[i]])<-rownames(rate.data)
    }
    do.call(rbind,rate.coef)->rate.coef
    colnames(rate.coef) <- c("slope","p-value")
    do.call(cbind,scalrat.data)->scalrat.data
    colnames(scalrat.data)<-colnames(rate.data)[1:(ncol(rate.data)-1)]
    data.frame(scalrat.data,age=rate.data$age)->scalrat.data

    rownames(rate.coef) <- names(rate.reg)<-colnames(rate.data)[1:(ncol(rate.data)-1)]

  }

  {##### Phenotypic Trend Real Multi #####
    phen.reg <-apply(phen.data[,1:(ncol(phen.data)- 1),drop=FALSE], 2, function(x) lm(range01(x) ~ phen.data[,ncol(phen.data)]))
    phen.coef <- lapply(phen.reg, function(x) coefficients(summary(x))[2,c(1,4)])
    # if(ncol(y)>1){
    #   if(is.null(colnames(y)))
    #     names(phen.coef) <- c(paste("y", seq(1:ncol(y)), sep = ""),"y.multiple") else
    #       names(phen.coef) <- c(colnames(y),"y.multiple")
    # }else names(phen.coef)<-"y"

    names(phen.coef) <- colnames(P)

    sapply(1:length(phen.coef),function(i){
      if(phen.coef[[i]][1]<0)
        (min(predict(phen.reg[[i]]))-mean(range01(phen.data[,i])))/sd(range01(phen.data[,i]))
      else (max(predict(phen.reg[[i]]))-mean(range01(phen.data[,i])))/sd(range01(phen.data[,i]))
    })->dev
    names(dev)<-names(phen.coef)
    do.call(rbind,phen.coef)->phen.coef
  }


  #### Nodes Real ####
  if (!is.null(node)) {
    rate.dataN <- rate.data
    phen.dataN <- phen.data
    rate.dataN$group <- rep("others", nrow(rate.dataN))
    phen.reg.sel<-phen.reg.age.sel<-phen.reg.y.sel<-list()
    rate.reg.y.sel<-rate.reg.age.sel<-list()
    for (j in 1:length(node)) {
      n <- node[j]
      sele <- getDescendants(t, n)
      sele[which(sele <= Ntip(t))]<-t$tip.label[sele[which(sele<=Ntip(t))]]
      rate.dataN[match(c(n,sele), rownames(rate.dataN)), ]$group <- paste("g",
                                                                          n, sep = "")
      rep("others",nrow(rate.dataN))->gg
      gg[match(c(n,sele), rownames(rate.dataN))]<-paste("g",n, sep = "")
      data.frame(rate.dataN,gg)->rate.dataN
      rate.data.sel <- rate.data[match(c(n,sele), rownames(rate.data)),,drop=FALSE]

      {#### Rate Trend Real Node Multi ####
        rate.reg.age <- rate.data.sel[,ncol(rate.data.sel),drop=FALSE]

        lapply(1:(ncol(rate.data.sel)-1),function(i){
          bet <- rate.data.sel[,i,drop=FALSE]
          age <- rate.data.sel[,ncol(rate.data.sel),drop=FALSE]
          abs(bet)->rts
          predict(lm(as.matrix(rts)~as.matrix(age)))->pp
          names(pp)<-rownames(rts)
          pp
        })->rate.reg.y

        # names(rate.reg.y) <- colnames(rate.data)[1:(ncol(rate.data)-1)]
        names(rate.reg.y) <- rownames(rate.coef)
      }
      rate.reg.y.sel[[j]] <- do.call(cbind,rate.reg.y)
      rate.reg.age.sel[[j]] <- rate.reg.age

      {##### Phenotypic Trend Real Node Multi #####
        PPsel <- phen.data[match(rownames(rate.data.sel), rownames(phen.data)),,drop=FALSE]
        phen.regN <- apply(PPsel[1:(ncol(PPsel)-1)],2, function(x)
          coefficients(summary(lm(range01(x) ~ PPsel[, ncol(PPsel)])))[2,c(1,4)])

        # if(ncol(y)>1){
        #   if(is.null(colnames(y)))
        #     colnames(phen.regN) <- c(paste("y", seq(1,dim(y)[2]), sep = ""),"y.multiple") else
        #       colnames(phen.regN) <- c(colnames(y),"y.multiple")
        # }else colnames(phen.regN)<-"y"

        phen.reg.y.sel[[j]] <- apply(PPsel[1:(ncol(PPsel)-1)], 2,function(x)
          predict(lm(x ~ PPsel[, ncol(PPsel)])))
        phen.reg.age.sel[[j]] <- PPsel$age
      }
      phen.reg.sel[[j]] <- phen.regN
    }

    names(phen.reg.sel)<-node

    #colnames(rate.dataN)[4:ncol(rate.dataN)]<-paste("group",seq(1,(ncol(rate.dataN)-3)),sep="")
    phen.dataN <- cbind(phen.dataN, rate.dataN[,(ncol(rate.data)+1):ncol(rate.dataN)])
    if (length(which(rate.dataN$group == "others")) < 3)
      rate.dataN <- rate.dataN[which(rate.dataN$group != "others"),]


    rate.dataRES<-rate.dataN[,1:(ncol(rate.dataN)-length(node))]

    { #### Node Comparison ####
      rate.NvsN<-rate.NvsO<-list()
      phen.NvsN<-phen.NvsN.emm<-phen.NvsO<-list()

      groupPP <- phen.dataN[which(phen.dataN$group != "others"), ]$group
      agePP <- phen.dataN$age

      for (i in 1:(ncol(rate.data)-1)) {
        bets <- rate.dataN[, i,drop=FALSE]
        yPP<-phen.dataN[,i,drop=FALSE]

        rate.comp<-list()
        phen.comp<-list()
        for (w in 1:(length(node))){
          data.frame(rate=unname(bets),age=rate.dataN$age,group=rate.dataN[,(w+ncol(rate.data)+1)])->emdat
          suppressMessages(mmeans<-as.data.frame(pairs(emmeans::emmeans(lm(abs(rate)~age+group,data=emdat),specs="group"))))
          # mtrends<-as.data.frame(pairs(emmeans::emtrends(lm(abs(rate)~age*group,data=emdat),specs="group",var="age")))
          # mtrends[match(mmeans[,1],mtrends[,1]),]->mtrends
          # data.frame(do.call(rbind,strsplit(as.character(mmeans[,1])," - ")),mmeans[,-c(1,3,4,5)],mtrends[,c(2,6)])->rate.comp[[w]]

          emmeans::emtrends(lm(abs(rate)~age*group,data=emdat),specs="group",var="age")->mtrends.test
          mtrends<-as.data.frame(pairs(mtrends.test))
          mtrends[match(mmeans[,1],mtrends[,1]),]->mtrends
          data.frame(do.call(rbind,strsplit(as.character(mmeans[,1])," - ")),mmeans[,-c(1,3,4,5)],
                     t(as.data.frame(mtrends.test)[,2,drop=FALSE]),mtrends[,6])->rate.comp[[w]]

          data.frame(yPP,age=phen.dataN$age,group=phen.dataN[,(w+ncol(rate.data)+1)])->PPemdat
          if((!is.null(x1))&isFALSE(x1.residuals)){ #### emmeans multiple ####
            data.frame(PPemdat,x1=x1[match(rownames(PPemdat),rownames(x1)),])->PPemdat
            suppressMessages(PPmeans<-as.data.frame(pairs(emmeans::emmeans(lm(range01(yPP[,1])~age+x1+group,data=PPemdat),specs="group"))))
          }else{
            suppressMessages(PPmeans<-as.data.frame(pairs(emmeans::emmeans(lm(range01(yPP[,1])~age+group,data=PPemdat),specs="group"))))
          }
          data.frame(do.call(rbind,strsplit(as.character(PPmeans[,1])," - ")),PPmeans[,-c(1,3,4,5)])->phen.comp[[w]]

        }
        do.call(rbind,rate.comp)->rate.NvsO[[i]]
        # colnames(rate.NvsO[[i]])<-c("group_1","group_2","emm.difference","p.emm","slope.difference","p.slope")
        colnames(rate.NvsO[[i]])<-c("group_1","group_2","emm.difference","p.emm","slope.node","slope.others","p.slope")

        do.call(rbind,phen.comp)->phen.comp
        colnames(phen.comp)<-c("group_1","group_2","mean","p.mean")
        sapply(strsplit(as.character(phen.comp[,1]),"g"),"[[",2)->PPnam
        phen.comp[,3:4]->phen.comp
        rownames(phen.comp)<-PPnam
        phen.comp[match(node,rownames(phen.comp)),]->phen.comp
        phen.comp->phen.NvsO[[i]]

        dat <- data.frame(bets=unname(bets), age=rate.dataN$age, group=rate.dataN$group)
        suppressMessages(mmeans<-as.data.frame(pairs(emmeans::emmeans(lm(abs(bets)~age+group,data=dat),specs="group"))))
        # mtrends<-as.data.frame(pairs(emmeans::emtrends(lm(abs(bets)~age*group,data=dat),specs="group",var="age")))
        # mtrends[match(mmeans[,1],mtrends[,1]),]->mtrends
        # data.frame(do.call(rbind,strsplit(as.character(mmeans[,1])," - ")),mmeans[,-c(1,3,4,5)],mtrends[,c(2,6)])->rate.NvsN[[i]]
        # colnames(rate.NvsN[[i]])<-c("group_1","group_2","emm.difference","p.emm","slope.difference","p.slope")
        # rate.NvsN[[i]][-which(rate.NvsN[[i]]$group_2=="others"),]->rate.NvsN[[i]]

        emmeans::emtrends(lm(abs(bets)~age*group,data=dat),specs="group",var="age")->mtrends.test
        as.data.frame(mtrends.test)[,1:2]->mtrends.slope
        mtrends<-as.data.frame(pairs(mtrends.test))
        mtrends[match(mmeans[,1],mtrends[,1]),]->mtrends
        data.frame(do.call(rbind,strsplit(as.character(mmeans[,1])," - ")),mmeans[,-c(1,3,4,5)])->mtrends.dat
        data.frame(mtrends.dat,mtrends.slope[match(mtrends.dat[,1],mtrends.slope[,1]),2],
                   mtrends.slope[match(mtrends.dat[,2],mtrends.slope[,1]),2],mtrends[,6])->rate.NvsN[[i]]
        colnames(rate.NvsN[[i]])<-c("group_1","group_2","emm.difference","p.emm","slope.group1","slope.group2","p.slope")
        rate.NvsN[[i]][which(rate.NvsN[[i]]$group_2!="others"),]->rate.NvsN[[i]]

        dat <- data.frame(yPP=unname(yPP), age=phen.dataN$age, group=phen.dataN$group)
        if((!is.null(x1))&isFALSE(x1.residuals)){ #### emmeans multiple ####
          data.frame(dat,x1=x1[match(rownames(dat),rownames(x1)),])->dat
          suppressMessages(PPpairs<-as.data.frame(pairs(emmeans::emmeans(lm(range01(yPP)~age+x1+group,data=dat),specs="group"))))
        }else{
          suppressMessages(PPpairs<-as.data.frame(pairs(emmeans::emmeans(lm(range01(yPP)~age+group,data=dat),specs="group"))))
        }

        data.frame(do.call(rbind,strsplit(as.character(PPpairs[,1])," - ")),
                   PPpairs[,-c(1,3,4,5)])->Pcomp.emm
        colnames(Pcomp.emm)<-c("group_1","group_2","mean","p.mean")
        Pcomp.emm[which(Pcomp.emm$group_2!="others"),]->Pcomp.emm
        Pcomp.emm->phen.NvsN.emm[[i]]

        if(length(node)>1){
          sapply(phen.reg.sel,function(w) w[1,i])->slope.tot
          names(slope.tot)<-paste("g",names(phen.reg.sel),sep="")

          combn(sort(unique(as.character(groupPP))),2)->pair

          # sapply(1:ncol(pair), function(jj){
          #   slope.tot[match(pair[1,jj],names(slope.tot))]-
          #     slope.tot[match(pair[2,jj],names(slope.tot))]
          # })->slope.diff
          # data.frame(t(pair),slope.diff)->phen.NvsN[[i]]
          # colnames(phen.NvsN[[i]])<-c(colnames(rate.NvsN[[i]])[1:2],"estimate")

          do.call(rbind,lapply(1:ncol(pair), function(jj){
            data.frame(slope.tot[match(pair[1,jj],names(slope.tot))],
                       slope.tot[match(pair[2,jj],names(slope.tot))])
          }))->slopes

          data.frame(t(pair),slopes)->phen.NvsN[[i]]
          colnames(phen.NvsN[[i]])<-c(colnames(rate.NvsN[[i]])[1:2],"slope.group_1","slope.group_2")

        }else{
          phen.NvsN<-NULL
          phen.NvsN.emm<-NULL
        }
      }

      if(!is.null(phen.NvsN.emm)) names(phen.NvsN.emm)<-colnames(phen.dataN)[1:ncol(rate.data)-1]

      sapply(strsplit(as.character(rate.NvsO[[1]][,1]),"g"),"[[",2)->grn
      lapply(rate.NvsO,function(x) x[match(node,grn),])->rate.NvsO

      lapply(1:length(node),function(u){
        # do.call(rbind,lapply(rate.NvsO,function(x) x[u,3:6]))->rcs
        do.call(rbind,lapply(rate.NvsO,function(x) x[u,3:7]))->rcs
        rownames(rcs)<-rownames(rate.coef)
        rcs
      })->rate.coef.sel

    }

    names(rate.coef.sel)<-node
    names(rate.NvsO) <- colnames(rate.dataN)[1:(ncol(rate.data)-1)]

    if(!is.null(phen.NvsN)){
      names(rate.NvsN)<-names(rate.NvsO)
      names(phen.NvsN.emm) <- names(phen.NvsN) <- names(phen.reg)
    }
  }


  #### Random ####
  RR$aces[1,,drop=FALSE]->a
  if(!is.null(x1)) x1[1:Ntip(t)]->y1
  {
    yyD <- list()
    yyT <- list()

    if(is.null(x1)){
      for (i in 1:ncol(y)) {
        yyD[[i]] <- suppressWarnings(replicate(nsim,fastBM(t, sig2 = 1, a = a[i], bounds=c(min(y[,i]),max(y[,i])))))
        yyT[[i]] <- replicate(nsim,fastBM(t, sig2 = 1, a = mean(y[,i]), bounds=c(min(y[,i]),max(y[,i]))))
      }
    }else{
      if(isFALSE(x1.residuals)){
        for (i in 1:ncol(y))
          yyD[[i]] <- suppressWarnings(replicate(nsim,fastBM(t, sig2 = 1, a = a[i], bounds=c(min(y[,i]),max(y[,i])))))
      }else{
        phen.data[which(rownames(phen.data)==(Ntip(t)+1)),]->ares
        phen.data[which(rownames(phen.data)%in%t$tip.label),]->yres
        for (i in 1:ncol(y))
          yyD[[i]] <- suppressWarnings(replicate(nsim,fastBM(t, sig2 = 1, a = ares[1,i], bounds=c(min(yres[,i]),max(yres[,i])))))
      }

      matrix(ncol=nsim,nrow=Ntip(t))->yy1T
      matrix(ncol=nsim,nrow=Nnode(t))->ace1T
      for(k in 1:nsim){
        phenb<-list()
        for(i in 1:ncol(y)) fastBM(t, sig2 = 1, a = mean(y[,i]), bounds=c(min(y[,i]),max(y[,i])),internal = TRUE)->phenb[[i]]
        fastBM(t,a=mean(y1),bounds=c(min(y1),max(y1)),internal = TRUE)->phen1b
        do.call(cbind,phenb)->phenb

        cor(cbind(y,y1))->cr
        cbind(phenb,phen1b)->xx1
        cr->m
        t(chol(m))->U
        U%*%t(xx1)->xx3
        xx3[-nrow(xx3),1:Ntip(t),drop=FALSE]->yyT[[k]]
        if(ncol(yyT[[k]])>nrow(yyT[[k]])) t(yyT[[k]])->yyT[[k]]
        xx3[nrow(xx3),1:Ntip(t)]->yy1T[,k]
        xx3[nrow(xx3),(Ntip(t)+1):(Ntip(t)+Nnode(t))]->ace1T[,k]
      }
      rownames(yy1T)<-rownames(y)
      rownames(ace1T)<-rownames(aceRR)
    }
  }

  ii=NULL
  res <- list()
  if(round((detectCores() * clus), 0)==0) cl<-makeCluster(1, setup_strategy = "sequential") else
    cl <- makeCluster(round((detectCores() * clus), 0), setup_strategy = "sequential")
  registerDoParallel(cl)
  res <- foreach(ii = 1:nsim, .packages = c("car","nlme", "ape",
                                            "phytools", "doParallel"
  )) %dopar% {
    # for(ii in 1:nsim){
    gc()
    #for(ii in 1:nsim){
    {
      yD <- lapply(yyD, function(x) x[,ii])
      yD <- do.call(cbind, yD)

      if(is.null(x1)){
        yT <- lapply(yyT, function(x) x[,ii])
        yT <- do.call(cbind, yT)
      }else{
        yT <-yyT[[ii]]
        y1T<-yy1T[,ii]
        aceT<-ace1T[,ii]
        cbind(L,y1T-aceT[1])->LX
        cbind(L1,aceT-mean(aceT))->L1X
      }

      betasT<- matrix(ncol=ncol(yT), nrow=Ntip(t)+Nnode(t))
      aceRRT<- matrix(ncol=ncol(yT), nrow=Nnode(t))
      betasD<- matrix(ncol=ncol(yD), nrow=Ntip(t)+Nnode(t))
      aceRRD<- matrix(ncol=ncol(yD), nrow=Nnode(t))
      for (s in 1:ncol(yT)){
        rootV<-a[s]
        lambda <- RR$lambda[s]

        if(!is.null(x1)){

          m.betasT <- (solve(t(LX) %*% LX + lambda * diag(ncol(LX))) %*%
                         t(LX)) %*% (as.matrix(yT[,s]) - rootV)
          aceRRT[,s] <- (L1X %*% m.betasT[c(1:Nnode(t),length(m.betasT)), ]) + rootV
          as.matrix(m.betasT[-length(m.betasT),])->betasT[,s]

        }else{
          m.betasT <- (solve(t(L) %*% L + lambda * diag(ncol(L))) %*%
                         t(L)) %*% (as.matrix(yT[,s]) - rootV)

          aceRRT[,s] <- (L1 %*% m.betasT[1:Nnode(t), ]) + rootV
          m.betasT->betasT[,s]
        }

        if(isTRUE(x1.residuals)) rootVD<-mean(range(yres[,s])) else rootVD<-mean(range(y[,s])) ### CONTROLLARE CON x1
        m.betasD <- (solve(t(L) %*% L + lambda * diag(ncol(L))) %*%
                       t(L)) %*% (as.matrix(yD[,s]) - rootVD)
        aceRRD[,s] <- (L1 %*% m.betasD[1:Nnode(t), ]) + rootVD
        m.betasD->betasD[,s]


      }
      rownames(betasT)<-rownames(betasD)<-rownames(RR$rates)
      rownames(aceRRT)<-rownames(aceRRD)<-rownames(RR$aces)

    }

    betasT->betasTreal
    #### Covariate ####
    if (!is.null(cov)) {
      cov[match(rownames(betas),names(cov))]->cov
      if (length(which(apply(betasT, 1, sum) == 0)) >
          0) {
        zeroes <- which(apply(betasT, 1, sum) == 0)
        R <- log(abs(betasT))
        R <- R[-zeroes, ]
        Y <- abs(cov)
        Y <- Y[-zeroes]
        resi <- residuals(lm(R ~ Y))
        factOut <- which(apply(betasT, 1, sum) != 0)
        betasT[factOut, ] <- resi
        betasT[zeroes, ] <- 0
      }else {
        R <- log(abs(betasT))
        Y <- abs(cov)
        resi <- residuals(lm(R ~ Y))
        betasT <- as.matrix(resi)
      }
      ratesT<-as.matrix(as.data.frame(apply(betasT, 1, function(x) sqrt(sum(x^2)))))
    }
    #### End of Covariate ####

    if (ncol(y)>1) {
      betasTreal->betasT
      ratesT<-as.matrix(as.data.frame(apply(betasT, 1, function(x) sqrt(sum(x^2)))))
      aceRRTmulti <- (L1 %*% ratesT[1:Nnode(t), ])
      yTmulti <- (L %*% ratesT)
      rate.data <- data.frame(betas = betasT[match(eds[, 1],rownames(betasT)), ],
                              rate = ratesT[match(eds[,1], rownames(ratesT)), ],
                              age = eds[, 2])
      colnames(rate.data)[1:ncol(yT)] <- paste("betas", seq(1,ncol(yT), 1), sep = "")

      ratesD<-as.matrix(as.data.frame(apply(betasD, 1, function(x) sqrt(sum(x^2)))))
      aceRRDmulti <- (L1 %*% ratesD[1:Nnode(t), ])
      yDmulti <- (L %*% ratesD)+aceRRDmulti[1,]

      nodesD<-cbind(aceRRD,aceRRDmulti)
      colnames(nodesD) <- c(paste("y", seq(1, dim(y)[2]),sep = ""),"y.multi")
      P <- rbind(nodesD, cbind(yD,yDmulti))
      if(isTRUE(x1.residuals)) apply(P,2,function(x) residuals(lm(x~x1)))->P

    }else {
      rate.data <- data.frame(rate = betasT[match(eds[, 1],rownames(betasT)), ],
                              age = eds[, 2])

      P <- rbind(aceRRD, yD)
      colnames(P)<-"y"
      if(isTRUE(x1.residuals)) as.matrix(residuals(lm(P~x1)))->P
    }
    rate.data[,ncol(rate.data)]+L[1,1]->rate.data[,ncol(rate.data)]
    phen.data <- data.frame(P[match(rownames(rate.data), rownames(P)), ],age=rate.data$age)

    {##### Rate Trend Random Multi #####
      rate.reg <- scalrat.data<-rate.coef <-list()
      ee<-array()
      for (i in 1:(ncol(rate.data)-1)) {
        bet <- rate.data[, i,drop=FALSE]
        age <- rate.data[, ncol(rate.data),drop=FALSE]
        abs(bet)->rts->rtsA
        log(range01(rts))->rts

        c(which(rts=="-Inf"),which(age=="-Inf"))->outs
        if(length(outs)>0){
          rts[-outs,,drop=FALSE]->rts
          age[-outs,,drop=FALSE]->age
        }
        age->ageC

        sd(range01(rtsA[ageC<0.5*max(ageC)]))/sd(range01(rtsA)[ageC>0.5*max(ageC)])->ee[i]

        if(!is.null(x1)){
          rts[-1,,drop=FALSE]->rts
          age[-1,,drop=FALSE]->age

          car::outlierTest(lm(as.matrix(rts)~as.matrix(age)))->ouT
          if(any(ouT$bonf.p<=0.05)){
            rts[-match(names(which(ouT$bonf.p<=0.05)),rownames(rts)),,drop=FALSE]->rts
            age[-match(names(which(ouT$bonf.p<=0.05)),rownames(age)),,drop=FALSE]->age
          }
        }else{
          lm(as.matrix(rts)~as.matrix(age))->bb

          residuals(bb)[order(residuals(bb),decreasing=TRUE)][1:(Ntip(t)/15)]->resout
          if((Ntip(t)+1)%in%names(resout)){
            rts[-1,,drop=FALSE]->rts
            age[-1,,drop=FALSE]->age
          }
        }
        lm(as.matrix(rts)~as.matrix(age))->regr.1
        rate.coef[[i]] <- coef(summary(regr.1))[2, c(1, 4)]
        rate.reg[[i]] <- regr.1

        scalrat.data[[i]]<-rts[match(rownames(rate.data),rownames(rts)),,drop=FALSE]
      }
      do.call(rbind,rate.coef)->rate.coef
      colnames(rate.coef) <- c("slope","p-value")
      do.call(cbind,scalrat.data)->scalrat.data
      colnames(scalrat.data)<-colnames(rate.data)[1:(ncol(rate.data)-1)]
      data.frame(scalrat.data,age=rate.data$age)->scalrat.data
      rownames(scalrat.data)<-rownames(rate.data)

      if(is.null(colnames(y))|ncol(y)==1) rownames(rate.coef) <- names(rate.reg)<-colnames(rate.data)[1:(ncol(rate.data)-1)] else
        rownames(rate.coef) <- names(rate.reg)<-c(colnames(y),"y.multiple")
    }


    { #### Phenotypic Trend Random Multi ####
      phen.reg <- apply(phen.data[1:(dim(phen.data)[2] - 1)], 2, function(x)
        summary(lm(range01(x) ~phen.data[, dim(phen.data)[2]])))
      phen.coef <- lapply(phen.reg, coefficients)
      if(ncol(y)>1){
        if(is.null(colnames(y)))
          names(phen.coef) <- c(paste("y", seq(1:ncol(y)), sep = ""),"y.multiple") else
            names(phen.coef) <- c(colnames(y),"y.multiple")
          cbind(yT,yTmulti)->yT
          cbind(yD,yDmulti)->yD
      }else names(phen.coef)<-"y"

    }


    #### Nodes Random ####
    if (!is.null(node)) {
      phen.dataN <- phen.data
      phen.dataN$group <- rep("NA", nrow(phen.dataN))
      phen.reg.sel <- list()
      for (j in 1:length(node)) {
        n <- node[j]
        sele <- getDescendants(t, n)
        sele[which(sele <= Ntip(t))]<-t$tip.label[sele[which(sele<=Ntip(t))]]
        phen.dataN[match(c(n,sele), rownames(phen.dataN)), ]$group <- paste("g",
                                                                            n, sep = "")
        {#### Phenotypic Trend Random Node Multi ####
          PPsel <- phen.data[match(c(n,sele), rownames(phen.data)),,drop=FALSE]
          phen.regN <- apply(PPsel[1:(ncol(PPsel)- 1)],
                             2, function(x) coefficients(summary(lm(range01(x) ~ PPsel[, ncol(PPsel)])))[2,c(1,4)])
        }
        phen.reg.sel[[j]] <- phen.regN
      }
      names(phen.reg.sel) <- node


      if(!is.null(phen.NvsN)){
        phen.dataN <- phen.dataN[-which(phen.dataN$group == "NA"),]

        {#### Node Comparison Random ####
          groupPP <- phen.dataN[, (ncol(phen.dataN))]
          agePP <- phen.dataN[, (ncol(phen.dataN)-1)]

          phen.NvsN.ran<-list()
          for (k in 1:(ncol(rate.data)-1)) {

            sapply(phen.reg.sel,function(w) w[1,k])->slope.tot
            names(slope.tot)<-paste("g",names(phen.reg.sel),sep="")


            combn(sort(unique(as.character(groupPP))),2)->pair
            sapply(1:ncol(pair), function(jj){
              slope.tot[match(pair[1,jj],names(slope.tot))]-
                slope.tot[match(pair[2,jj],names(slope.tot))]
            })->slope.diff

            data.frame(t(pair),slope.diff)->phen.NvsN.ran[[k]]
            colnames(phen.NvsN.ran[[k]])<-c("group_1","group_2","estimate")

          }
          names(phen.NvsN.ran)<-colnames(phen.dataN)[1:(dim(phen.dataN)[2] - 2)]
        }
      }else{
        phen.NvsN.ran<-NULL
      }
      res[[ii]] <- list(phen.coef, rate.coef,
                        phen.data, rate.data,yT,yD,ee,scalrat.data,
                        phen.reg.sel,phen.NvsN.ran)
    }else {
      res[[ii]] <- list(phen.coef, rate.coef,
                        phen.data, rate.data,yT,yD,ee,scalrat.data)
    }
  }
  stopCluster(cl)
  #### End of Random ####

  {#### p Rate Trend Multi####
    p.rate<-spread<-array()
    rate.coefR<-list()
    scalrat.CI <- CIabsolute <- list()
    if(ncol(y)>1) cbind(y,y.multi)->yTot else y->yTot
    for (i in 1:(ncol(rate.data)-1)){
      scalrat.CI[[i]]<-rate.coefR <- do.call(rbind, lapply(lapply(res,"[[", 2), function(x) x[i, 1]))
      # quantile(rate.coefR,c(0.025,0.975))
      yRT<- lapply(lapply(res,"[[",5), function(x) x[, i])

      rank(c(rate.coef[i, 1],rate.coefR[-1]))[1]/nsim->pp->p.rate[i]
      sapply(lapply(res,"[[",7),"[[",i)->ee

      if(i<=ncol(y)){
        ifelse(rate.coef[i, 1]>0,{
          if(pp<=0.025) dir.rate<-"increase slower" else
            if(pp>=0.975) dir.rate<-"increase faster" else
              dir.rate<-"nosig"
        },{
          if(pp>=0.975) dir.rate<-"decrease slower" else
            if(pp<=0.025) dir.rate<-"decrease faster" else
              dir.rate<-"nosig"
        })

        if(dir.rate!="nosig"){
          if(ncol(y)>1) yprint<-paste(" for",colnames(phen.data)[i]) else yprint<-NULL
          print(paste("Absolute evolutionary rates ",
                      dir.rate," than BM expectation",yprint,sep=""))
        }

      }

      (diff(range(yTot[,i]))/diff(range(yTot[,i][which(diag(vcv(t))<diff(range(diag(vcv(t)))))])))*
        mean(sapply(yRT,function(x)(diff(range(x))/diff(range(x[which(diag(vcv(t))<diff(range(diag(vcv(t)))))])))))->spread[i]

      as.matrix(sapply(lapply(res,"[[", 4),function(k) unname(coef(lm(abs(k[,i])~k$age))[2])))->CIabsolute[[i]]

    }
    names(spread)<-names(p.rate) <- rownames(rate.coef)
    p.rate <- data.frame(slope = rate.coef[, 1],p.real = rate.coef[, 2],
                         p.random = p.rate,spread=spread)
  }

  { #### p Phenotypic Trend Multi and pdf ####
    CIphenotype<-list()
    p.phen <- array()
    for (i in 1:(ncol(phen.data)-1)) {
      CIphenotype[[i]]<-sapply(lapply(lapply(res,"[[", 1), "[[", i), function(x) x[2,1])
      p.phen[i] <- rank(c(phen.coef[i,1],CIphenotype[[i]][-1]))[1]/nsim

      if(i<=ncol(y)){
        ifelse(phen.coef[i,1]>0,{
          if(p.phen[i]<=0.025) dir.phen<-"increases less" else
            if(p.phen[i]>=0.975) dir.phen<-"increases more" else
              dir.phen<-"nosig"
        },{
          if(p.phen[i]<=0.025) dir.phen<-"decreases more" else
            if(p.phen[i]>=0.975) dir.phen<-"decreases less"else
              dir.phen<-"nosig"
        })

        if(dir.phen!="nosig"){
          if(ncol(y)>1) yprint<-paste(" for",colnames(phen.data)[i]) else yprint<-NULL
          print(paste("Phenotypic regression slope ",dir.phen ," than BM expectation",yprint,sep=""))
        }
      }
    }

    p.phen <- cbind(phen.coef, p.phen,dev)
    colnames(p.phen) <- c("slope", "p.real",
                          "p.random","dev")

  }

  if (!is.null(node)) {
    p.phen.sel <- list()
    for (u in 1:length(node)) {
      {#### p Phenotypic Trend Node Multi ####
        p.sele <- list()
        for (i in 1:(ncol(phen.data)-1)) {
          slopeR <- sapply(lapply(lapply(res,"[[", 9), "[[", u),function(a) a[1,i])
          p.sel <- rank(c(phen.reg.sel[[u]][1,i],slopeR[-1]))[1]/nsim

          p.sele[[i]] <- c(slope = phen.reg.sel[[u]][1,i],
                           p.slope = p.sel,
                           emm.difference=phen.NvsO[[i]][match(names(phen.reg.sel)[u],rownames(phen.NvsO[[i]])),1],
                           p.emm=phen.NvsO[[i]][match(names(phen.reg.sel)[u],rownames(phen.NvsO[[i]])),2])

        }
        names(p.sele) <- colnames(phen.reg.sel[[u]])
        p.phen.sel[[u]]  <- do.call(rbind, p.sele)
      }
    }
    names(p.phen.sel) <- names(phen.reg.sel)


    p.rate.sel<-rate.coef.sel
    p.rate.comp <- rate.NvsN
    if(!is.null(phen.NvsN)){
      {#### p Phenotypic Slope Comparison Multi  ####
        p.phen.comp<-list()
        for (k in 1:(ncol(phen.data)-1)) {
          p.nodeR<-mean.diff<-array()
          for(i in 1:nrow(phen.NvsN[[k]])){
            rank(c(phen.NvsN[[k]][i,3]-phen.NvsN[[k]][i,4],sapply(lapply(lapply(res,"[[",10),"[[",k),function(x) x[i,3])[-1]))[1]/nsim->p.nodeR[i]
            if((phen.NvsN[[k]][i,3]-phen.NvsN[[k]][i,4])>0&p.nodeR[i]<=0.05|(phen.NvsN[[k]][i,3]-phen.NvsN[[k]][i,4])<0&p.nodeR[i]>=0.95) p.nodeR[i]<-0.5
            if(as.character(interaction(phen.NvsN[[k]][i,c(1,2)]))==as.character(interaction(phen.NvsN.emm[[k]][i,c(1,2)])))
              phen.NvsN.emm[[k]][i,3]->mean.diff[i] else
                -1*phen.NvsN.emm[[k]][i,3]->mean.diff[i]
          }
          # data.frame(phen.NvsN[[k]][,c(1,2)],slope.difference=phen.NvsN[[k]][,3],p.slope=1-p.nodeR,
          #            emm.difference=mean.diff,p.emm=phen.NvsN.emm[[k]][,4])->p.phen.comp[[k]]
          data.frame(phen.NvsN[[k]][,c(1,2)],slope.group_1=phen.NvsN[[k]][,3],slope.group_2=phen.NvsN[[k]][,4],p.slope=p.nodeR,
                     emm.difference=mean.diff,p.emm=phen.NvsN.emm[[k]][,4])->p.phen.comp[[k]]
        }
        names(p.phen.comp)<-names(phen.reg)
      }
    }else{
      p.phen.comp<-NULL
    }

    data.frame(scalrat.data,group=rate.dataRES[match(rownames(scalrat.data),rownames(rate.dataRES)),]$group)->scalrat.data
  }

  CInts <- list(CIphenotype, CIabsolute,scalrat.CI)
  names(CInts) <- c("phenotype", "absolute_rate","rescaled_rate")
  class(CInts)<-"RRphyloList"

  list(phen.dataRES,rate.dataRES,scalrat.data)->tabs
  names(tabs) <- c("phenotypeVStime", "absrateVStime","rescaledrateVStime")
  class(tabs)<-"RRphyloList"

  if (!is.null(node)) {
    {
      for (i in 1:length(p.rate.sel)) {
        for(k in 1:ncol(y)){
          ifelse(p.phen.sel[[i]][k,1]>0,{
            if(p.phen.sel[[i]][k,2]<=0.025) phenN.dir<-" increases less" else
              if(p.phen.sel[[i]][k,2]>=0.975) phenN.dir<-" increases more" else
                phenN.dir<-"nosig"
          },{
            if(p.phen.sel[[i]][k,2]<=0.025) phenN.dir<-" decreases more" else
              if(p.phen.sel[[i]][k,2]>=0.975) phenN.dir<-" decreases less" else
                phenN.dir<-"nosig"
          })

          if(phenN.dir!="nosig"){
            if(ncol(y)>1) yprint<-paste(" for",rownames(p.phen.sel[[i]])[k]) else yprint<-NULL
            print(paste("Phenotypic regression slope through node ",names(p.rate.sel)[i],
                        phenN.dir," than BM expectation",yprint,sep=""))
          }


          if(p.rate.sel[[i]][k,3]>p.rate.sel[[i]][k,4]) rateN.dir<-"higher" else rateN.dir<-"lower"
          if(p.rate.sel[[i]][k,5]<=0.05){
            if(ncol(y)>1) yprint<-paste(" for",rownames(p.phen.sel[[i]])[k]) else yprint<-NULL
            print(paste("Absolute evolutionary rates regression slope",yprint,
                        " through node ",names(p.rate.sel)[i],
                        " is ",rateN.dir," than for the rest of the tree",sep=""))
          }
        }
      }

      if(!is.null(p.phen.comp)){
        for (j in 1:ncol(y)) {
          for (i in 1:nrow(p.phen.comp[[j]])) {
            if ((p.phen.comp[[j]][i,3]-p.phen.comp[[j]][i,4])>0&p.phen.comp[[j]][i,5]>=0.95) phen.comp.dir<-"higher" else
              if ((p.phen.comp[[j]][i,3]-p.phen.comp[[j]][i,4])<0&p.phen.comp[[j]][i,5]<=0.05) phen.comp.dir<-"lower" else
                phen.comp.dir<-NULL
              if(!is.null(phen.comp.dir)){
                gsub("g","",p.phen.comp[[j]][i,1:2])->nprint
                if(ncol(y)>1) yprint<-paste(" for",names(p.phen.comp)[j]) else yprint<-NULL
                print(paste("Phenotypic regression through node ",nprint[1],
                            " is ",phen.comp.dir," than regression through node ",nprint[2],
                            yprint,sep=""))
              }

              if (p.rate.comp[[j]][i, 7] <= 0.05){
                if (p.rate.comp[[j]][i, 5]<p.rate.comp[[j]][i, 6]) rate.comp.dir<-"lower" else rate.comp.dir<-"higher"
                gsub("g","",p.rate.comp[[j]][i,1:2])->nprint
                if(ncol(y)>1) yprint<-paste(" for",names(p.rate.comp)[j]) else yprint<-NULL
                print(paste("Absolute evolutionary rates regression slope through node ",nprint[1],
                            " is ",rate.comp.dir,"than regression through node ",nprint[2],yprint,sep=""))
              }
          }
        }
      }
    }

    if(is.null(p.phen.comp)){
      node.comp.res<-NULL
    }else{
      if(length(p.phen.comp)==1) p.phen.comp[[1]]->p.phen.comp
      if(length(p.rate.comp)==1) p.rate.comp[[1]]->p.rate.comp
      node.comp.res <- list(p.phen.comp, p.rate.comp)
      names(node.comp.res) <- c("phenotype", "rate")
    }

    res <- list(tabs, p.phen, p.rate,
                p.phen.sel, p.rate.sel,
                node.comp.res, CInts)
    names(res) <- c("trend.data", "phenotypic.regression", "rate.regression",
                    "node.phenotypic.regression", "node.rate.regression",
                    "group.comparison", "ConfInts")



  }else {
    res <- list(tabs, p.phen, p.rate,
                CInts)
    names(res) <- c("trend.data", "phenotypic.regression", "rate.regression",
                    "ConfInts")
  }

  if(length(which(sapply(res,function(x) is.null(x))))>0)
    res<-res[-which(sapply(res,function(x) is.null(x)))]
  return(res)
}