## podar arbol
opar <- par
library(ape)
hm288 <- read.tree("/home/m/Dropbox/Phylogenies & data/hummer_corrected_names_288sp.tre")
hm288$tip.label <- gsub("_", " ", hm288$tip.label)
s <- hm288$tip.label
hm288$tip.label <- paste(toupper(substring(s, 1, 1)), substring(s, 2), sep = "")
drs <- list.dirs(path = "/home/m/Desktop/grabaciones", full.names = F)
drs <- drs[-1]
drs <- drs[drs %in% hm288$tip.label]
drs <- drs[drs != "Heliothryx barroti"]
selsp <- c(drs, grep("leucotis", hm288$tip.label, value = T), c("Hylocharis eliciae", "Hylocharis cyanus"))
hm19 <- drop.tip(hm288, tip = hm288$tip.label[!hm288$tip.label %in% selsp])
hm19$tip.label[hm19$tip.label == "Hylocharis leucotis"] <- "Basilinna leucotis"
plot(hm19)
# est.rep <- sample(c("lek", "no.lek"), 19, replace = T)
#
# complx <- rep(0, length(est.rep))
#
# complx[est.rep == "lek"] <- rnorm(n = length(which(est.rep == "lek")), mean = 10, sd = 3)
#
# complx[est.rep != "lek"] <- rnorm(n = length(which(est.rep != "lek")), mean = 7, sd = 3)
# library(phytools)
# phylANOVA(tree = hm19, x = est.rep, y = complx)
#
# complx2 <- rnorm(n = 19, mean = 10, sd = 2)
#
# phylANOVA(tree = hm19, x = est.rep, y = complx2)
#
ANOVA filogenetica
library(phytools)
## Loading required package: maps
dat <- read.csv("/home/m/Dropbox/Sharing/Evo canto colibries OET/Datos complejidad canto.csv")
##
dat2 <- aggregate(dat[, 3:5], by = list(dat$especie), mean)
names(dat2)[1] <- "especie"
dat2 <- merge(dat2, dat[!duplicated(dat$especie),c(1, 6)], by = "especie", all = F)
dat2 <- dat2[match(hm19$tip.label, dat2$especie),]
datsd <- aggregate(dat[, 3:5], by = list(dat$especie), sd)
names(datsd)[1] <- "especie"
datsd <- datsd[match(hm19$tip.label, datsd$especie),]
names(datsd)[2:4] <- paste0(names(datsd)[2:4], "sd")
rownames(dat2) <- dat2$especie
dat3 <- data.frame(dat2, datsd[,-1])
estrategia <- dat2$estrategia
names(estrategia) <- dat2$especie
n.categorias <- dat2$n.categorias
names(n.categorias) <- dat2$especie
phylANOVA(tree = hm19, x = estrategia, y = n.categorias, nsim = 1000)
## $F
## [1] 9.102936
##
## $Pf
## [1] 0.001
##
## $T
## Lek No.lek
## Lek 0.000000 -3.017107
## No.lek 3.017107 0.000000
##
## $method
## [1] "holm"
##
## $Pt
## Lek No.lek
## Lek 1.000 0.001
## No.lek 0.001 1.000
# boxplot(n.categorias ~ estrategia, ylab = "Categorias de elementos", xlab = "Estrategia reproductiva", col = "pink")
dat4 <- aggregate(dat3[,c(2:4,6:8)],by = list(dat3$estrategia), mean)
names(dat4)[1] <- "estrategia"
library(ggplot2)
ggplot(dat4, aes(x= estrategia, y=n.categorias)) +
geom_line() +
geom_errorbar(width=.1, aes(ymin=n.categorias-n.categoriassd, ymax=n.categorias+n.categoriassd)) + ylab("Número de categorias") +
geom_point(shape=21, size=3, fill="white")
## geom_path: Each group consists of only one observation. Do you need to
## adjust the group aesthetic?
n.elementos <- dat2$n.elementos
names(n.elementos) <- dat2$especie
phylANOVA(tree = hm19, x = estrategia, y = n.elementos, nsim = 1000)
## $F
## [1] 0.3857531
##
## $Pf
## [1] 0.353
##
## $T
## Lek No.lek
## Lek 0.0000000 -0.6210902
## No.lek 0.6210902 0.0000000
##
## $method
## [1] "holm"
##
## $Pt
## Lek No.lek
## Lek 1.000 0.353
## No.lek 0.353 1.000
# boxplot(n.elementos ~ estrategia, ylab = "Número de elementos", xlab = "Estrategia reproductiva", col = "pink")
ggplot(dat4, aes(x= estrategia, y=n.elementos)) +
geom_line() +
geom_errorbar(width=.1, aes(ymin=n.elementos-n.elementossd, ymax=n.elementos+n.elementossd)) + ylab("Número de elementos") +
geom_point(shape=21, size=3, fill="white")
## geom_path: Each group consists of only one observation. Do you need to
## adjust the group aesthetic?
duracion <- dat2$duracion
names(duracion) <- dat2$especie
phylANOVA(tree = hm19, x = estrategia, y = duracion, nsim = 1000)
## $F
## [1] 0.2074986
##
## $Pf
## [1] 0.49
##
## $T
## Lek No.lek
## Lek 0.0000000 -0.4555202
## No.lek 0.4555202 0.0000000
##
## $method
## [1] "holm"
##
## $Pt
## Lek No.lek
## Lek 1.00 0.49
## No.lek 0.49 1.00
ggplot(dat4, aes(x= estrategia, y=duracion)) +
geom_line() +
geom_errorbar(width=.1, aes(ymin=duracion-duracionsd, ymax=duracion+duracionsd)) + ylab("Duración (s)") +
geom_point(shape=21, size=3, fill="white")
## geom_path: Each group consists of only one observation. Do you need to
## adjust the group aesthetic?
Señal filogenetica
phylosig(tree = hm19, method = "lambda", x = n.categorias, nsim = 1000, test = T)
## $lambda
## [1] 6.933942e-05
##
## $logL
## [1] -15.86082
##
## $logL0
## [1] -15.86067
##
## $P
## [1] 1
phylosig(tree = hm19, method = "lambda", x = n.elementos, nsim = 1000, test = T)
## $lambda
## [1] 1.048675
##
## $logL
## [1] -52.3953
##
## $logL0
## [1] -55.73984
##
## $P
## [1] 0.009700441
phylosig(tree = hm19, method = "lambda", x = duracion, nsim = 1000, test = T)
## $lambda
## [1] 6.933942e-05
##
## $logL
## [1] -39.56552
##
## $logL0
## [1] -39.56538
##
## $P
## [1] 1
par(mar = c(5, 4, 10, 2) + 0.1)
b<-as.matrix(dat2[,2:4])
colnames(b) <- c("n.cat", "n.elem", "dur")
phylo.heatmap(tree = hm19,X = b,lwd=3,standardize=TRUE, fsize = c(1, 0.7, 1))
Polimorphoespacios
par(mar = c(5, 4, 10, 2) + 0.1)
colors <- setNames(c("blue", "green"),c("Lek", "No.lek"))
dotTree(hm19, estrategia, colors = colors)
hm19a <- hm19
sp <- dat2$especie
sp <- as.character(sp)
sp[estrategia == "Lek"] <- paste0(sp[estrategia == "Lek"], "*")
hm19a$tip.label <- sp
n.elementos2 <- dat2$n.elementos
names(n.elementos2) <- sp
hm19a <- chronoMPL(hm19a)
par(mar = c(5, 4, 4, 2) + 0.1)
phenogram(hm19a, n.elementos2, spread.labels=TRUE, colors = "blue", ylab = "Número de elementos")
n.categorias2 <- dat2$n.categorias
names(n.categorias2) <- sp
phenogram(hm19a, n.categorias2, spread.labels=TRUE, colors = "blue", ylab = "Número de categorias", xlab = "Tiempo")
duracion <- dat2$duracion
names(duracion) <- sp
phenogram(hm19a, duracion, spread.labels=TRUE, colors = "blue", ylab = "Duración", xlab = "Tiempo")
