rm(list=ls())
# Load necessary libraries
library(ape)

# Add taxa I before (A,B)
tree_string <- "((A,(B,C))G,(F,(D,E))H);"
tree <- read.tree(text = tree_string)

# Define edge colors and widths
edge_colors <- c(rep("blue", tree$Nnode), rep("red", tree$Nnode))
edge_widths <- c(rep(5, tree$Nnode), rep(5, tree$Nnode))

plot(tree, edge.color = edge_colors, edge.width = edge_widths, cex=1.5, label.offset = 0.5)
title(expression(paste("Two AR(1) (", phi[1], ", ", phi[2], ") vs. One AR(1) Rate (", phi, ")")))

#plot(tree, edge.color = edge_colors, edge.width = edge_widths, cex=1.5,  label.offset = 0.5)
#title("Two ARs vs. One AR Rate")
text(x = 0.05, y = 5.35, adj = c(0, 0.5), labels = expression(paste(sigma[des], " = ", phi[1], sigma[anc], " + ", epsilon[des - anc])), col = "red")
text(x = 1.25, y = 3.35, adj = c(0, 0.5), labels = expression(paste(sigma[des], " = ", phi, sigma[anc], " + ", epsilon[des - anc])))
text(x = 0.05, y = 1.05, adj = c(0, 0.5), labels = expression(paste(sigma[des], " = ", phi[2], sigma[anc],  " + ", epsilon[des - anc])), col = "blue")