rm(list=ls()) install.packages("ape") install.packages("geiger") install.packages("phytools") install.packages("caper") install.packages("graphics") library(ape) library(geiger) library(phytools) library(caper) library(graphics) library(phytools) X<-read.table("rawdata_Baeckens_Whiting_ProcB_2020.txt", header=T,sep="\t") head(X) tree.hac<-read.nexus("pyron_tree.nex") #reads tree tree<- drop.tip(tree.hac,tip=tree.hac$tip.label[!tree.hac$tip.label %in% X$Species]) tree ###PHYLOGENETIC SIGNAL compdata<-comparative.data(phy=tree,data=X,names.col=Species,vcv=FALSE,vcv.dim=3,na.omit=FALSE) phylo.d(compdata, binvar=FP) phylo.d(compdata, binvar=SG) ###FITPAGEL TEST X<-read.table("rawdata_Baeckens_Whiting_ProcB_2020.txt", header=T,sep="\t",row.names=1) head(X) x<-setNames(X$FP,rownames(X)) y<-setNames(X$SG,rownames(X)) head(x) head(y) fit.xy<-fitPagel(tree,x,y) fit.xy fit.x<-fitPagel(tree,x,y,dep.var="x") fit.x fit.y<-fitPagel(tree,x,y,dep.var="y") fit.y aic<-setNames(c(fit.xy$independent.AIC, fit.x$dependent.AIC, fit.y$dependent.AIC, fit.xy$dependent.AIC), c("independent","dependent x", "dependent y","dependent x&y")) aic aic.w(aic) ####MuSSE library(diversitree) tree2<-force.ultrametric(tree,method="extend") st<-setNames(X$state,rownames(X)) lik.m <- make.musse(tree2, st, 4) argnames(lik.m) #UNCONSTRAINTED lik.base <- constrain(lik.m, q14 ~ 0, q41 ~ 0, q23 ~ 0, q32~ 0) p <- starting.point.musse(tree2, 4) fit.base <- find.mle(lik.base, p[argnames(lik.base)]) fit.base #CONSTRAINED lik.lambda <- constrain(lik.m, q14 ~ 0, q41 ~ 0, q23 ~ 0, q32~ 0, q12 ~ q34) fit.lambda <- find.mle(lik.lambda, p[argnames(lik.lambda)]) fit.lambda anova(fit.base, lambda = fit.lambda) #inverse prediction: sociality promotes the evolution of epidermal glands? lik.lambda2 <- constrain(lik.m, q14 ~ 0, q41 ~ 0, q23 ~ 0, q32~ 0, q13 ~ q24) fit.lambda2 <- find.mle(lik.lambda2, p[argnames(lik.lambda2)]) fit.lambda2 anova(fit.base, lambda = fit.lambda2) ###ANC X<-read.table("rawdata_Baeckens_Whiting_ProcB_2020.txt", header=T,sep="\t",row.names=1) head(X) x<-setNames(X$FP,rownames(X)) y<-setNames(X$SG,rownames(X)) head(x) head(y) mx<-make.simmap(tree,x,nsim=100) my<-make.simmap(tree,y,nsim=100) overlap<-sapply(mx,function(x,y) sapply(y,map.overlap,x),y=my) dim(overlap) mean(overlap) hist(overlap,xlim=c(0,1),col="grey") lines(c(mean(overlap),mean(overlap)),c(0,par()$usr[4]),lty="dashed",col="red",lwd=2) abline(v=0.5,lty="dashed") m<-mean(overlap) head(m) head(overlap) overlap t.test(overlap,mu=0.5)