bm.multiphylo <- function(trees, matrix, subsample = NULL, stype = NULL, delta, criterion, Barplot = FALSE) {
 
  ## Testa se os pacotes requeridos estão instalados (código modificado de https://stackoverflow.com/questions/9341635/check-for-installed-packages-before-running-install-packages)
  for (i in c("ape", "geiger")) { # Cria um ciclo que itera sobre o vetor com os nomes dos pacotes 
    if (i %in% installed.packages()) { # Procura cada um dos pacotes na lista de pacotes instalados no computador 
      cat (i, "is already installed in this computer\n") # Se encontrar o pacote procurado, imprime na tela uma mensagen
    } else {
        stop (i, "is not installed\n") # Se não encontra o pacote, para e imprime na tela uma mensagen 
      }
  }
 
  ## Carrega os pacotes requeridos
  library(ape)
  library(geiger)
 
  ## Leia os dados morfológicos (o objeto `matrix`) e extrai o número de dimensões
  if (ncol(matrix) > 2) { # Testa se o objeto `matrix` tem mais de dois colunas
    stop ("The object `matrix` must contain only two columns, or you must
       to specify the column to be analyzed\n") # Se o objeto `matrix` tem mais de dois colunas, para e imprime uma mensagen na tela
  } else {
      chtaxa <- as.character(matrix[ ,1]) # Cria um vetor com os nomes dos táxons inclusos no objeto `matrix`  
    }
 
  ## Leia o arquivo multiphylo contendo as árvores filogenéticas (o objeto `trees`)
  if (class(trees) == "multiPhylo") { # Testa se o objeto `trees` é de clase multiPhylo
    trtaxa <- trees$TipLabel$tip.label # Cria um vetor com os nomes dos táxons inclusos no objeto `trees`
  } else {
      stop ("The object `trees` must be of class multiphylo, with a sample of trees
            in parenthetical notation\n") # Para a função se o objeto `trees` não é de clase multiPhylo 
    }
 
  ## Testa se as árvores fologenéticas no objeto `trees` estão enraizadas
  if (min(is.rooted(trees)) == 0) { # Testa se as árvores filogenéticas estão enraizadas
    stop ("The phylogenetic trees must be rooted\n") # Para a função se pelo menos uma das árvores filogenéticas não está enraizada
  }
 
  ## Amostras de árvores filogenéticas
  if (!(is.null(subsample))) { # Testa se o usuário deu um valor para o arg `subsample`
    if (subsample >= 2 & subsample < length(trees)) { # Testa se o valor do argumento subsample está entre 2 e o número total de árvores
      stype <- match.arg(arg = stype, choices = c("random", "first", "last")) # Estabelece os valores possíveis do argumento `stype` 
      if (stype == "random") {
        trees <- sample(x = trees, size = subsample) # Se stype = "random" o subsample é ao acaso
      }
      if (stype == "first") {
        trees <- trees[1:subsample] # Se stype = "first" o subsample extrai as n primeiras árvores
      } 
      if (stype == "last") {
        trees <- trees[(length(trees) - subsample):(length(trees) -1)] # Se stype = "last" o subsample extrai as n últimas árvores
      }  
 
    } else {
      stop ("subsample must be numeric between 2 and the total of trees in the object trees -1 \n") # Se o valor do argumento `subsample` está errado imprime na tela uma mensagem
    }  
  }
 
  ## Compara os nomes dos táxons que estão na matriz morfológica (`matrix`) e nas árvores filogenéticas (`trees`)
  if (sum(sort(chtaxa) != sort(trtaxa)) > 0){ # Testa se os táxons nos objetos `matrix` e `trees` são iguais, para isso primeiro ordena os nomes alfabeticamente
 
    cat ("The following taxa in `matrix` were not found in the object `trees`.
         They will be removed for the analysis:\n",
         chtaxa[chtaxa %in% trtaxa == FALSE], "\n") # Imprime na tela os nomes dos táxons que estão no objeto `matrix` mas não no objeto `trees`
 
    cat ("The following taxa in `trees` were not found in the object `matrix`.
         They will be removed for the analysis:\n",
         trtaxa[trtaxa %in% chtaxa == FALSE], "\n") # Imprime na tela os nomes dos táxons que estão no objeto `trees` mas não no objeto `matrix`
 
    ## Remove os táxons incompativeis da `matrix`
    chtaxa.remove <- chtaxa[chtaxa %in% trtaxa == FALSE] # Cria um vetor com os nomes dos táxons que estão no objeto `matrix` mas não no objeto `trees`
 
    for (i in 1:length(chtaxa.remove)){ # Cria um contador desde 1 até o número de elementos no vetor `chtaxa.remove` 
      matrix <- matrix[-which(matrix[, 1] == chtaxa.remove[i]), ] # Remove do objeto `matrix` cada um dos táxons que não estão no objeto `trees`
    }
 
    ## Remove os táxons incompativeis do objeto `trees` (Código modificado de http://blog.phytools.org/)
    trees <- lapply (X = trees, FUN = drop.tip, tip = c(trtaxa[trtaxa %in% chtaxa == FALSE])) # Aplica a função `drop.tip` sobre cada uma das árvores do objeto `trees`
    class(trees) <- "multiPhylo" # Muda a clase do novo objeto `trees` de lista para multiphylo
 
    cat ("trees:", length(row.names(summary(trees))), "phylogenetic trees with",
         length(trees[[1]]$tip.label), "taxa\n") # Imprime na tela o número de árvores no objeto `trees` e o novo número de táxons que contem cada uma das árvores
  } else {
      ## Sumário dos dados que serão analisados
      # Árvores filogenéticas
      cat ("trees:", length(row.names(summary(trees))), "phylogenetic trees with",
           length(trees$TipLabel$tip.label), "taxa\n") # Imprime na tela o número de árvores no objeto `trees` e o número de táxons que contem cada uma das árvores
    } 
 
  ## Sumário dos dados que serão analisados
  # Matriz de dados morfológicos
  matrix[ ,2] <- as.factor(matrix[ ,2]) # Muda a clase do carater analisado
  cat ("matrix: 1 character and", dim(matrix)[[1]], "coded taxa\n") # 
  ch <- setNames(object = matrix[ , 2], nm = matrix[ , 1]) # Cria um objeto contendo só os estados do carater
  cat ("The character has the following states:\n", levels(ch), "\n") # Extrai a lista de estados do carater e os imprime na tela
 
  ## Use a função `fitDiscrete` para testar o melhor modelo de evolução de caracteres para cada uma das árvores do objecto `trees`
  models <- c("ER", "ARD", "SYM") # Cria um objeto com os nomes dos modelos de evolução que vão ser testados
  fitmodels <- data.frame(Tree = rep(1:length(trees), each = 3), # Cria um data frame para inserir o resultado que sai do for
                          Model = rep(c("ER", "ARD", "SYM"), times = length(trees)),
                          AIC = rep(NA, times = (length(trees))*3),
                          AICc = rep(NA, times = (length(trees))*3))
 
  iter <- 1 # Cria um objeto que aumenta em 1 com cada ciclo, permite inserir resultados nas colunas do data frame
  for (i in 1:length(trees)) { # Cria um contador de 1 até o número de árvores no objeto `trees`
    cat("Running", i, "trees out of", length(trees),"\n") # Imprime na tela uma mensagem com o número de ciclos por vez
    for(j in models) { # Cria um contador que itera sobre o vetor dos modelos
      tree <- trees[[i]] # Seleciona uma árvore do objeto `trees`
      testfit <- fitDiscrete(phy = tree, dat = ch, model = j) # Aplica a função `fitDiscrete` sobre a árvore selecionada
      fitmodels[iter, "AIC"] <- testfit$opt$aic # Extrai o valor de AIC do objeto testfit e o salva no data frame
      fitmodels[iter, "AICc"] <- testfit$opt$aicc # Extrai o valor de AICc do objeto testfit e o salva no data frame
      iter <- iter + 1 # Aumenta em 1 o valor do objeto `iter`. Permite inserir resultados na próxima linha do data frame
    }
  }
 
  ## Calcula a diferença (delta) entre os modelos de evolução
  minAIC <- aggregate(fitmodels$AIC ~ fitmodels$Tree, FUN = min) # Extrai o mínimo valor de AIC para cada árvore
  minAIC <- rep(minAIC$`fitmodels$AIC`, each = 3) # Cria um vetor repitindo mínimo valor de AIC para cada árvore
 
  minAICc <- aggregate(fitmodels$AICc ~ fitmodels$Tree, FUN = min) # Extrai o mínimo valor de AICc para cada árvore
  minAICc <- rep(minAICc$`fitmodels$AICc`, each = 3) # Cria um vetor repitindo mínimo valor de AICc para cada árvore
 
  fitmodels$delta_AIC <- fitmodels$AIC - minAIC # Calcula a diferença (delta) entre o menor valor de AIC e os restantes; salva o resultado numa nova coluna do data frame `firmodels`
  fitmodels$delta_AICc <- fitmodels$AICc - minAICc # Calcula a diferença (delta) entre o menor valor de AICc e os restantes; salva o resultado numa nova coluna do data frame `firmodels`
 
  ## Seleciona o melhor modelo usando os criterios AIC e AICc
  bestModel_AIC <- fitmodels[which(fitmodels$delta_AIC == 0), c("Tree", "Model")] # Extrai o nome do modelo com o valor mais baixo de AIC (o melhor modelo) para cada árvore
  bestModel_AICc <- fitmodels[(which(fitmodels$delta_AICc == 0)), c("Tree", "Model")] # Extrai o nome do modelo com o valor mais baixo de AIC (o melhor modelo) para cada árvore
 
  ## Calcula os "Akaike weights" para cada um dos modelos usando o valor de AICc (Uma pior versão da função aic.w de `phytools`)
  fitmodels$Tree <- as.factor(fitmodels$Tree) # Muda a clase da coluna Tree do data frame `fitmodels`
  b <- rep(NA, times = length(levels(fitmodels$Tree))*3) # Cria um vetor para salvar o resultado do for
  iter = 1 # # Cria um objeto que aumenta em 1 com cada ciclo, permite inserir resultados no vetor `b`
  for (i in levels(fitmodels$Tree)) { # Cria um contador que itera sobre o número de árvores
    for (j in models) { # Cria um contador que itera sobre os modelos
      a <- fitmodels[which(fitmodels$Tree == i & fitmodels$Model == j), "AICc"] # Seleciona o valor de AICc para cada árvore e cada modelo, e o salva no objeto `a`
      b[iter] <- (exp(1))^(- a/2) # Calcula uma parte da expresão matemática que calcula os pesos dos modelos (ver documentação)
      iter = iter + 1 # Aumenta em 1 o valor do objeto `iter`. Permite inserir resultados no vetor `b`
    }
  }
 
  fitmodels$b <- b # Cria uma nova coluna no data frame `fitmodels` para salvar o conteudo do vetor `b`
  sumation <- aggregate(x = fitmodels$b, by = list(Tree = fitmodels$Tree), FUN = sum) # Soma os valores de `b` para cada árvore
  sumation <- rep(sumation$x, each = 3) # Repite os resultados da soma anterior para criar um vetor do mesmo comprimento que as colunas no data frame `fitmodels`
  fitmodels$AICcw <- b / sumation # Calcula os pesos dos modelos (ver documentação) e salva-los numa nova coluna do data frame `fitmodels` 
 
  ## Seleciona o melhor modelo usando o criterio AICcw
  max_AICcw <- aggregate(x = fitmodels$AICcw, by = list(Tree = fitmodels$Tree), FUN = max) # Extrai o maior valor de AICcw (o melhor modelo) para cada árvore
  bestModel_AICcw <- fitmodels[which(fitmodels$AICcw %in% max_AICcw$x), 1:2] # Extrai o nome do modelo com o maior valor AICcw (o melhor modelo) para cada árvore  
 
  ## Data frame com o melhor modelo selecionado para cada árvore
  bestModels <- data.frame(Tree = levels(fitmodels$Tree), AIC = bestModel_AIC$Model, # Cria um data frame com os nomes dos melhores modelos para cada árvore de acordo a cada um dos criterios
                           AICc = bestModel_AICc$Model, AICcw = bestModel_AICcw$Model)
 
  ## Re-organiza o data frame fitmodels
  fitmodels <- data.frame(Tree = fitmodels$Tree, Model = fitmodels$Model, AIC = fitmodels$AIC, # Reordena as colunas do data frame fitmodels
                          AICc = fitmodels$AICc, AICcw = fitmodels$AICcw, delta_AIC = fitmodels$delta_AIC,
                          delta_AICc = fitmodels$delta_AICc)
 
  ## Compara os valores de AIC e AICc entre modelos e seleciona um modelo de acordo com o valor do arg `delta`
  if (class(delta) == "numeric") { # Testa se a clase do objeto `delta` é numeric
    deltaAIC <- fitmodels[which(fitmodels$delta_AIC <= delta), c(1,2,6)] # Seleciona o modelo que tem um valor de delta menor ou igual ao valor de delta_AIC inserido pelo usuário
    deltaAICc <- fitmodels[which(fitmodels$delta_AICc <= delta), c(1,2,7)] # Seleciona o modelo que tem um valor de delta menor ou igual ao valor de delta_AICc inserido pelo usuário
 
    ## Troca o nome de cada modelo por o número relativo de parámetros de cada um deles para AIC
    deltaAIC$param <- rep(NA, times = length(deltaAIC$Model)) # Cria uma nova coluna no data frame `deltaAIC` para salvar os resultados do for
    models <- c("ER", "SYM", "ARD") # # Cria um objeto com os nomes dos modelos de evolução ordenados pela sua complexidade
    for (i in 1:3) { # Cria um contador de 1 até o número de modelos
      deltaAIC$param[which(deltaAIC$Model == models[i])] = i # Escreve 1, 2 ou 3 na nova coluna do `deltaAIC` para cada um dos modelos de acordo com sua complexidade
    }
 
    min_AICparam <- aggregate(x = deltaAIC$param, by = list(Tree = deltaAIC$Tree), FUN = min) # Extrai o menor valor de parámetros (o modelo mais ssimples) para cada árvore
    selectModel_AIC <- deltaAIC[which(deltaAIC$param %in% min_AICparam$x), c("Tree", "Model")] # Extrai o nome do modelo mais simples de acorco com o AIC para cada árvore  
 
    ## Troca o nome de cada modelo por o número relativo de parámetros de cada um deles para AICc
    deltaAICc$param <- rep(NA, times = length(deltaAICc$Model)) # Cria uma nova coluna no data frame `deltaAIC` para salvar os resultados do for
    for (i in 1:3) { # Cria um contador de 1 até o número de modelos
      deltaAICc$param[which(deltaAICc$Model == models[i])] = i # Escreve 1, 2 ou 3 na nova coluna do `deltaAIC` para cada um dos modelos de acordo com sua complexidade
    }
 
    min_AICcparam <- aggregate(x = deltaAICc$param, by = list(Tree = deltaAICc$Tree), FUN = min) # Extrai o menor valor de parámetros (o modelo mais ssimples) para cada árvore
    selectModel_AICc <- deltaAICc[which(deltaAICc$param %in% min_AICcparam$x), c("Tree", "Model")] # Extrai o nome do modelo mais simples de acorco com o AICc para cada árvore
 
  } else {
      stop ("arg `delta` must be numeric") # Se o valor do argumento `delta` inserido pelo usuário não é numérico, para a função e imprime uma mensagem na tela
    }
 
  ## Sumárico gráfico
  if (Barplot == TRUE) { # Testa se o usuário quer gerar uma saída gráfica
    bestModels$Tree <- as.factor(bestModels$Tree) # Confere que a coluna Tree seja de clase factor
    b <- c(as.character(bestModels$AIC), as.character(bestModels$AICc), as.character(bestModels$AICcw)) # Cria um vetor com o melhor modelo para cada árvore e cada criterio
    c <- data.frame(criterion = rep(c("AIC", "AICc", "AICcw"), each = length(levels(bestModels$Tree))), 
                    models = b)
    freqModels <- table(c) # Cria uma tabela para contar as veces que cada modelo foi o melhor modelo ajustado para a amostra de árvores
 
    x11() # Abre um dispositivo gráfico
    par(mfrow = c(2,1), mar = c(5, 6, 4, 4))
    color.names = c("black","grey40","grey80") # Estabelece as cores de cada categoría para o gráfico
    barplot(t(freqModels), beside = TRUE, ylim = c(0, max(freqModels)+1), xlab = "Akaike criteria for model selection",
            ylab = "Model Frequency", col = color.names, axis.lty = "solid", main = "Best Model using Akaike criteria") # Faz um gráfico de barplot para resumir a contagem dos modelos
    legend("topright", legend = rownames(t(freqModels)), fill = color.names, title = "Models", cex = 0.5) # Insere uma legenda
 
    criterion <- match.arg(arg = criterion, choices = c("AIC", "AICc", "AICcw")) # Estabelece os valores possíveis do argumento `criterion`
    if (criterion == "AIC") {
      freqModels2 <- table(selectModel_AIC$Model)
    }
 
    if (criterion == "AICc") {
      freqModels2 <- table(selectModel_AICc$Model)
    }
 
    if (criterion == "AICcw") {
      freqModels2 <- table(bestModel_AICcw$Model)
    }
 
    barplot(freqModels2, beside = TRUE, ylim = c(0, max(freqModels2)+1), main = "Selected Model" , xlab = "Evaluated model for character evolution",
            ylab = "Model Frequency", col = color.names, axis.lty = "solid") # Faz um gráfico de barplot para resumir a contagem dos modelos
  }  
 
  ## Saída
  final_AIC <- list(matrix = matrix, trees = trees, fitModels = fitmodels, bestModels = bestModels,
                    selectModel = selectModel_AIC) # Cria  uma lista com os elementos que serão retornados pela função, caso o usuário use `criterion` = AIC
 
  final_AICc <- list(matrix = matrix, trees = trees, fitModels = fitmodels, bestModels = bestModels,
                     selectModel = selectModel_AICc) # Cria  uma lista com os elementos que serão retornados pela função, caso o usuário use `criterion` = AICc
 
  final_AICcw <- list(matrix = matrix, trees = trees, fitmodels = fitmodels, bestModels = bestModels) # Cria  uma lista com os elementos que serão retornados pela função, caso o usuário use `criterion` = AICcw
 
  if (criterion == "AIC") return(final_AIC) # Para a função e retorna uma das lista já ditas
  if (criterion == "AICc") return(final_AICc) # Para a função e retorna uma das lista já ditas
  if (criterion == "AICcw") return(final_AICcw) # Para a função e retorna uma das lista já ditas
}

Função function.r

Documentação documentation.txt

Dados morfológicos morphodata.csv

Árvores filogenéticas (mudar extensão a .tre) multiphylo.txt