library(devtools) install_github("willgearty/pcmtools")

library(pcmtools)

library(OUwie) data(tworegime) ou.results <- list() ou.results[[1]] <- OUwie(tree,trait,model=c("BM1")) ou.results[[2]] <- OUwie(tree,trait,model=c("BMS"), root.station = FALSE) ou.results[[3]] <- OUwie(tree,trait,model=c("OUM")) ou.results[[4]] <- OUwie(tree,trait,model=c("OUMV")) #Both regimes have same parameters for BM1 model. Both regimes have different parameters for other models. regime.mat <- data.frame(BM1 = c(1, 1), BMS = c(1,2), OUM = c(1,2), OUMV = c(1,2), row.names = c(1,2)) #Extract and map parameters to regimes ou.parameters <- OUwieParSumm(ou.results, regime.mat) #Summarize model fit ou.aic <- OUwieAICSumm(ou.parameters) #Model average parameters ou.avg <- OUwieModelAvg(ou.parameters)

library(phytools) #Simulate a mapped tree set.seed(4) Q <- matrix(c(-2,1,1,1,-2,1,1,1,-2),3,3) rownames(Q) <- colnames(Q) <- letters[1:3] tree <- sim.history(pbtree(n=100,scale=1),Q) cols <- setNames(c("blue","red","green","orange"),letters[1:4]) #Plot the mapping plot(tree, cols, ftype="i", fsize=0.7) #Split state c to state d within subclade tree2 <- mergeMappedStates2(tree, "c", "d", 173) #Plot new mapping plot(tree2, cols, ftype="i", fsize=0.7)

library(phytools) #Simulate tree tree <- pbtree(n=70,scale=1) #Simulate discrete trait Q <- matrix(c(-1,1,1,-1),2,2) rownames(Q) <- colnames(Q)<-c(0,1) x1 <- sim.history(tree,Q)$states #Generate stochastic maps and density map mtrees <- make.simmap(tree,x1,nsim=100) map1 <- densityMap(mtrees) #Simulate continuous trait x2 <- fastBM(tree,sig2=0.1) #Generate cont map map2 <- contMap(tree, x2) #Extract posterior densities pp_data <- extractSimmapDensity(map1, map2) #Plot to see how traits have evolved with respect to one another plot(pp_data$map1, pp_data$map2)

library(phytools) #Simulate tree tree <- pbtree(n=70,scale=1) #Simulate discrete trait Q <- matrix(c(-1,1,1,-1),2,2) rownames(Q) <- colnames(Q)<-c(0,1) tree <- sim.history(tree,Q) #Simulate continuous trait x2 <- fastBM(tree,sig2=0.1) #Calculate stats through time tstt <- traitStatsThroughTime(tree, x2) coords <- tstt$xy stats <- tstt$stats #Plot phenogram ys <- c(coords$y1, coords$y2) plot(NULL, xlim=c(0,1), ylim=c(min(ys), max(ys)), ylab="trait", xlab="time") segments(coords$x1, coords$y1, coords$x2, coords$y2) #Plot stats plot(stats$x, stats$mean)