Reconstrução reticulada de áreas ancestrais (Dispersal-Extinction Cladogenesis model)
José Ricardo Inacio Ribeiro
Isso é um R Markdown Notebook depositado no Bioforum. Quando executar o código, os resultados aparecerão abaixo dele!
Agora, vamos implementar um modelo mais complexo. Vamos inserir eventos anagenéticos de dispersão e extinção no modelo.
## Biogeografia paramétrica:
library(ape)
library(phangorn)
library(phytools)## Carregando pacotes exigidos: mapslibrary(geiger)
#library(gtools)
## ler árvore de uma string:
texto.arv <- ('((x1:0.5, x2:0.5):0);')
texto.arv <- ('((x1:10, x2:10):10, x3:20);')
arv1 <- read.tree(text = texto.arv)
plot(arv1, no.margin = TRUE, edge.width = 2)
# vamos ver os nós:
#plotaremos com os nós terminais e nós internos x4 e x5:
plotTree(arv1, offset = 1, type="phylogram")
labelnodes(1:(Ntip(arv1) + arv1$Nnode),
1:(Ntip(arv1) + arv1$Nnode),
interactive = F,cex = 0.8)
Implementei uma função alterada que leva em consideração os processos de anagênese, chamada ZECA.biogeo_DEC, para três ranges A, B e AB.
# n.range = número de letras usadas
# n.areas = número máximo de ranges nos nós
ZECA.biogeo.DEC <- function(ranges, n.range, n.areas, arvore, x.D, x.E){
arv <- arvore
# áreas A, B e C e a matriz de transição Q (forma analítica):
#x.D <- taxaD
#x.E <- taxaE
Parea <- 1 / (2^n.range - 1) # priori de cada área (pi)
priori <- matrix(c(Parea, 0, 0,
0, Parea, 0,
0, 0, Parea), nc = 3, byrow = T)
#Cladogênese (sem o parâmetro j):
cladogen <- matrix(c(1/3, 0, 1/18,
0, 1/3, 1/18,
1/18, 1/18, 1/18), byrow = T, ncol = 3)
colnames(cladogen) <- rownames(cladogen) <- c('A', 'B', 'AB')
arranjo <- t(combn(c(LETTERS[1:n.range]), n.areas))
for(j in 1:dim(arranjo)[1]){
for(i in 1:n.range){
if (all(ranges == LETTERS[i])){
letras <- which(LETTERS[1:n.range] != LETTERS[i])
ranges[length(ranges) + 1] <- paste0(arranjo[j,], collapse = '')
ranges[length(ranges) + 1] <- LETTERS[letras]
}
}
if (any(ranges != LETTERS[1])){
ranges[length(ranges) + 1] <- paste0(arranjo[j,], collapse = '')
}
if (all(ranges == paste0(arranjo[j,], collapse = ''))){
ranges <- c(ranges, LETTERS[1:n.range])
}
}
tras <- matrix(c(NA, 0, x.D, 0, NA, x.D, x.E, x.E, NA),
byrow = TRUE, ncol = 3)
#trans.sup <- upper.tri(trans, diag = FALSE) * x.D
#trans.inf <- lower.tri(trans, diag = FALSE) * x.E
#trans <- trans.sup + trans.inf
diag(tras) <- -rowSums(tras, na.rm = T)
names(ranges) <- c(arv$tip.label, ranges[(length(arv$tip.label) + 1):length(ranges)])
if (!is.factor(ranges)) {
ranges <- factor(ranges)
}
x <- reorder(x = as.factor(ranges), X = c(1:length(ranges)))
nll <- nlevels(x)
cores <- palette()
cor <- cores[as.integer(x)]
cols <- setNames(cor, x)
ARE <- to.matrix(ranges, x)
ARE <- ARE[-c(dim(ARE)[1]:length(arv$tip.label) + 1),]
##################################
# algoritmo para cálculo de verossimilhanças:
lvls <- levels(x)
x <- as.integer(x)
rate <- tras #coeficientes da matriz de transição
#likelihood
nb.tip <- length(arv$tip.label)
if(any(arvore$edge.length == 0)){
nb.node <- arv$Nnode - 1
} else{
nb.node <- arv$Nnode
}
liks <- matrix(0, nb.tip + nb.node, nll) # matriz que irá possuir as probabilidades marginais condicionais
veros <- matrix(0, nb.tip + nb.node, nll) # matriz que irá possuir as verossimilhanças fracionadas
# priors <- matrix(0, nb.tip + nb.node, nll) # suposta matriz de prioris de áreas em um nó
TIPS <- 1:nb.tip
liks[cbind(TIPS, x[1:nb.tip])] <- 1
# priors[cbind(TIPS, x[1:nb.tip])] <- 1
#caso não haja estado no terminal
if (anyNA(ranges)) {
liks[which(is.na(x)), ] <- 1
# priors[which(is.na(x)), ] <- 1
}
arv <- reorder(arv, "postorder") #reordenando os levels
e1 <- arv$edge[, 1]
e2 <- arv$edge[, 2]
EL <- arv$edge.length
# e1 #nós internos
# e2 #nós terminais e raiz
# EL #ramos
# optamos pela decomposição espectral (uma forma analítica):
comp <- numeric(nb.tip + nb.node)
Q <- tras #matriz de transição de ranges
decompo <- eigen(Q)
lambda <- decompo$values
GAMMA <- decompo$vectors
invGAMMA <- solve(GAMMA)
#para a área
for (i in seq(from = 1, by = 2, length.out = nb.node)) {
j <- i + 1L
if(any(arvore$edge.length == 0)){
anc <- e1[i] - 1
} else {
anc <- e1[i]
}
des1 <- e2[i]
des2 <- e2[j]
#prioris = estimando-se as probabilidades da evolução dos
#ranges ancestral - descendente (anagênese):
dQi <- GAMMA %*% diag(exp(lambda * EL[i])) %*% invGAMMA
# p.anai <- priori %*% dQi * liks[des1,] #prioris para a dispersão e extinção
dQj <- GAMMA %*% diag(exp(lambda * EL[j])) %*% invGAMMA
# p.anaj <- priori %*% dQj * liks[des2,] #prioris para a dispersão e extinção
#probabilidades de cada range
# prioriT_i <- priori %*% cladogen %*% dQi
prioriT <- priori %*% cladogen * 1
# prioriT_j <- priori %*% cladogen %*% dQj
# Verossimilhanças condicionais dos processos de anagênese e da cladogênese
v.l <- liks[des1,] %*% dQi %*% prioriT
# v.l <- liks[des1,] %*% p.anai %*% prioriT_i
v.r <- liks[des2,] %*% dQj %*% prioriT
v <- v.l * v.r # verossimilhança do cenário: L(nodes)
veros[anc,] <- v
comp[anc] <- sum(v)
liks[anc, ] <- v / comp[anc]
}
#UFA! As reconstruções marginais rapidinhas!
x <- reorder(x = as.factor(ranges), X = c(1:length(ranges)))
colnames(liks) <- levels(x)
colnames(veros) <- levels(x)
#forma matricial:
vero.total <- veros[anc,] %*% prioriT
vero.total <- sum(vero.total) # verossimilhança do modelo
#log(vero.total) #em log na base e
#################################################
# RECONSTRUÇÃO:
plotTree(arv, type = "phylogram", fsize = 0.7, ftype = "i",lwd = 1, offset = 2)
if(any(arvore$edge.length == 0)){
nodelabels(node = 1:(arv$Nnode - 1) + Ntip(arv) + 1,
pie = liks[1:(arv$Nnode - 1) + Ntip(arv),],
piecol = c('black', '#61D04F', '#DF536B'), cex = 0.8)
} else {
nodelabels(node = 1:arv$Nnode + Ntip(arv),
pie = liks[1:arv$Nnode + Ntip(arv),],
piecol = unique(cols), cex = 0.8)
}
ARE <- to.matrix(ranges, levels(x))
legend(x = "bottomleft", legend = unique(x),
pch = 19, col = unique(cols), bty = 'o',
pt.cex = 1.5)
tiplabels(pie = ARE[arv$tip.label,], piecol = unique(cols),
cex = 0.5)
#######################
probMar <- liks[-TIPS,]
if(all(arvore$edge.length != 0)){
rownames(probMar) <- c(anc, anc + 1)
}
verosCond <- veros[-TIPS,]
if(all(arvore$edge.length != 0)){
rownames(verosCond) <- c(anc, anc + 1)
}
return(list(prob.marginal = probMar, vero_fracionada = verosCond, log_vero = log(vero.total)))
}Vamos comparar alguns modelos já conhecidos nossos?
\(\color{blue}{\text{Condição geral: a taxa total é baixa. Modelo 1: $D < E$}}\)
# the overall rate is low
# D < E
mod1 <- ZECA.biogeo.DEC(ranges = c('A', 'B', 'B'), n.range = 2, n.areas = 2, arvore = arv1,
x.D = 0.001, x.E = 0.009)
mod1 #k = 2 parâmetros livres (2^n - 2, com n = 2)## $prob.marginal
## A B AB
## 4 0.003966072 0.89675581 0.09927812
## 5 0.060829314 0.06082931 0.87834137
##
## $vero_fracionada
## A B AB
## 4 1.369997e-05 3.097656e-03 0.0003429355
## 5 2.374992e-05 2.374992e-05 0.0003429355
##
## $log_vero
## [1] -7.769615\(\color{blue}{\text{modelo 2: $D = E$}}\)
# D = E
mod2 <- ZECA.biogeo.DEC(ranges = c('A', 'B', 'B'), n.range = 2, n.areas = 2, arvore = arv1,
x.D = 0.005, x.E = 0.005)
mod2 #k = 2 parâmetros livres## $prob.marginal
## A B AB
## 4 0.01818511 0.8942056 0.08760933
## 5 0.19056280 0.1905628 0.61887440
##
## $vero_fracionada
## A B AB
## 4 7.118329e-05 0.0035002534 0.0003429355
## 5 1.055962e-04 0.0001055962 0.0003429355
##
## $log_vero
## [1] -7.637533\(\color{blue}{\text{modelo 3: $D > E$}}\)
# D > E
mod3 <- ZECA.biogeo.DEC(ranges = c('A', 'B', 'B'), n.range = 2, n.areas = 2, arvore = arv1,
x.D = 0.009, x.E = 0.001)
mod3 #k = 2 parâmetros livres## $prob.marginal
## A B AB
## 4 0.02912014 0.8823741 0.0885058
## 5 0.24719462 0.2471946 0.5056108
##
## $vero_fracionada
## A B AB
## 4 0.0001128325 0.0034189558 0.0003429355
## 5 0.0001676622 0.0001676622 0.0003429355
##
## $log_vero
## [1] -7.648253# the overall rate is high
# D < E
mod11 <- ZECA.biogeo.DEC(ranges = c('A', 'B', 'B'), n.range = 2, n.areas = 2, arvore = arv1,
x.D = 0.01, x.E = 0.09)
mod11 #k = 2 parâmetros livres (2^n - 2, com n = 2)## $prob.marginal
## A B AB
## 4 0.07170083 0.8705725 0.05772671
## 5 0.34397646 0.3439765 0.31204708
##
## $vero_fracionada
## A B AB
## 4 0.0004259512 0.0051717863 0.0003429355
## 5 0.0003780255 0.0003780255 0.0003429355
##
## $log_vero
## [1] -7.20255# D = E
mod22 <- ZECA.biogeo.DEC(ranges = c('A', 'B', 'B'), n.range = 2, n.areas = 2, arvore = arv1,
x.D = 0.05, x.E = 0.05)
mod22 #k = 2 parâmetros livres## $prob.marginal
## A B AB
## 4 0.2453614 0.6742048 0.08043386
## 5 0.4253696 0.4253696 0.14926087
##
## $vero_fracionada
## A B AB
## 4 0.0010461158 0.0028745202 0.0003429355
## 5 0.0009773113 0.0009773113 0.0003429355
##
## $log_vero
## [1] -7.547774# D > E
mod33 <- ZECA.biogeo.DEC(ranges = c('A', 'B', 'B'), n.range = 2, n.areas = 2, arvore = arv1,
x.D = 0.09, x.E = 0.01)
mod33 #k = 2 parâmetros livres## $prob.marginal
## A B AB
## 4 0.2873506 0.5590167 0.1536327
## 5 0.4058044 0.4058044 0.1883913
##
## $vero_fracionada
## A B AB
## 4 0.0006414175 0.0012478245 0.0003429355
## 5 0.0007387005 0.0007387005 0.0003429355
##
## $log_vero
## [1] -8.239737\(\color{blue}{\text{CONCLUSÃO:}}\) Se os ranges dos terminais forem, em sua maioria, ranges de áreas únicas, os dados serão geralmente muito bem ajustados a esse modelo se as taxa de extinção ou de dispersão forem muito altas! **Lembre-se: quanto maior for a taxa de extinção, menor será o tempo em que ranges amplos persistirão na evolução do grupo de estudo e menor será a probabilidade dos terminais estarem em um range amplo. A melhor forma de se evitar o range nulo é diminuindo a taxa de extinção!
“Model behavior is data-dependent, there is not ‘one’ behavior!”
(MASSANA et al., 2015)
Avaliando-se os modelos com log-likelihood ratio test:
Lembrando que o número de parâmetros livres, k, é igual a: \(k = 2\)
############# Avaliação dos modelos ##############
# Vamos comparar modelos que são aninhados (likelihood-ratio tests).
modelos<-c("mod1", "mod2", "mod3", "mod22", "mod33")
graus <- c(2, 2, 2, 2, 2, 2)
LR.test<-matrix(NA, 5, 3, dimnames = list(paste(modelos, "-vs-mod11", sep = ""),
c("LR", "df", "P(X2)")))
allfits<-list(mod1$log_vero, mod2$log_vero, mod3$log_vero,
mod22$log_vero, mod33$log_vero, mod11$log_vero)
for(i in 2:length(allfits) - 1){
LR.test[i, 1] <- - 2 * (allfits[[i]] - allfits[[6]])
LR.test[i, 2] <- graus[i] - graus[6]
LR.test[i, 3]<-pchisq(LR.test[i, 1], df = LR.test[i, 2],
lower.tail = FALSE)
}
LR.test## LR df P(X2)
## mod1-vs-mod11 1.1341304 0 0
## mod2-vs-mod11 0.8699670 0 0
## mod3-vs-mod11 0.8914070 0 0
## mod22-vs-mod11 0.6904484 0 0
## mod33-vs-mod11 2.0743740 0 0-Isso significa que a variação na taxa de extinção ou na taxa de dispersão pode melhorar ou piorar o ajuste dos dados!
Avaliando-se os modelos com \(AIC\), \(AIC_c\) e \(BIC\), estatísticas que não partem da premissa que os modelos devam ser aninhados…
# Valor de AIC (e os pesos de AKAIKE) do nossos modelos ajustados e
# não aninhados:
modelos <- c('mod1', 'mod2', 'mod3', 'mod11', 'mod22', 'mod33')
k <- c(2, 2, 2, 2, 2, 2)
n <- 1 # número de caracteres (Butler and King, 2004)
stats.test<-matrix(NA, 6, 12, dimnames = list(paste(modelos), c("LR", "df", "P(X2)",
"AIC", "d.AIC", "wAIC", "AICc","d.AICc", "wAICc", "BIC","d.BIC",
"wBIC")))
AIC.val1 <- -2 * mod1$log_vero + (2 * k[1]) #AIC
AIC.val2 <- -2 * mod2$log_vero + (2 * k[2]) #AIC
AIC.val3 <- -2 * mod3$log_vero + (2 * k[3]) #AIC
AIC.val11 <- -2 * mod11$log_vero + (2 * k[4]) #AIC
AIC.val22 <- -2 * mod22$log_vero + (2 * k[5]) #AIC
AIC.val33 <- -2 * mod33$log_vero + (2 * k[6]) #AIC
AIC.valores <- c(AIC.val1, AIC.val2, AIC.val3, AIC.val11, AIC.val22, AIC.val33)
round(AIC.valores, 3)## [1] 19.539 19.275 19.297 18.405 19.096 20.479# AICc:
AICc.val1 <- AIC.val1 * (n / (n - k[1] - 1))
AICc.val2 <- AIC.val2 * (n / (n - k[2] - 1))
AICc.val3 <- AIC.val3 * (n / (n - k[3] - 1))
AICc.val11 <- AIC.val11 * (n / (n - k[4] - 1))
AICc.val22 <- AIC.val22 * (n / (n - k[5] - 1))
AICc.val33 <- AIC.val33 * (n / (n - k[6] - 1))
AICc.valores <- c(AICc.val1, AICc.val2, AICc.val3, AICc.val11, AICc.val22, AICc.val33)
round(AICc.valores, 3)## [1] -9.770 -9.638 -9.648 -9.203 -9.548 -10.240# Valores de BIC (Bayesian Information Criterion):
BIC.val1 <- -2 * mod1$log_vero + (k[1] * log(n))
BIC.val2 <- -2 * mod2$log_vero + (k[2] * log(n))
BIC.val3 <- -2 * mod3$log_vero + (k[3] * log(n))
BIC.val11 <- -2 * mod11$log_vero + (k[4] * log(n))
BIC.val22 <- -2 * mod22$log_vero + (k[5] * log(n))
BIC.val33 <- -2 * mod33$log_vero + (k[6] * log(n))
BIC.valores <- c(BIC.val1, BIC.val2, BIC.val3, BIC.val11, BIC.val22, BIC.val33)
round(BIC.valores, 3)## [1] 15.539 15.275 15.297 14.405 15.096 16.479w.aic <- 0
delta.aic <- 0
w.aicc <- 0
delta.aicc <- 0
w.bic <- 0
delta.bic <- 0
for(i in 1:length(k)){
delta.aic[i] <- AIC.valores[i] - min(AIC.valores)
delta.aicc[i] <- AICc.valores[i] - min(AICc.valores)
delta.bic[i] <- BIC.valores[i] - min(BIC.valores)
}
w.aic<- (exp(- delta.aic / 2)) / (exp( - delta.aic[1] / 2) +
exp( - delta.aic[2] / 2) + exp( - delta.aic[3] / 2))
w.aicc<- (exp(- delta.aicc / 2)) / (exp( - delta.aicc[1] / 2) +
exp( - delta.aicc[2] / 2) + exp( - delta.aicc[3] / 2))
w.bic<- (exp(- delta.bic / 2)) / (exp( - delta.bic[1] / 2) +
exp( - delta.bic[2] / 2) + exp( - delta.bic[3] / 2))
stats.test[,4] <- round(AIC.valores, 3)
stats.test[,5] <- round(delta.aic, 3)
stats.test[,6] <- round(w.aic, 3)
stats.test[,1] <- c(round(as.numeric(LR.test[1:3,1]), 3), '-', round(as.numeric(LR.test[4:5,1]), 3))
stats.test[,2] <- c(round(as.numeric(LR.test[1:3, 2]), 3), '-', round(as.numeric(LR.test[4:5, 2]), 3))
stats.test[,3] <- c(round(as.numeric(LR.test[1:3, 3]), 3), '-', round(as.numeric(LR.test[4:5, 3]), 3))
stats.test[,7] <- round(AICc.valores, 3)
stats.test[,8] <- round(delta.aicc, 3)
stats.test[,9] <- round(w.aicc, 3)
stats.test[,10] <- round(BIC.valores, 3)
stats.test[,11] <- round(delta.bic, 3)
stats.test[,12] <- round(w.bic, 3)
stats.test <- as.data.frame(stats.test) #resultados
stats.test