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