library(tidyselect)
library(PhylogeneticEM)
library(viridis)

forelimb_lsr <- as.matrix(select(data, ends_with("_lsr_fore"))) #select log shape ratio traits for forelimbs


######### forelimb Standard PCA ######### 
forelimb_lsr_pca <- princomp(forelimb_lsr, cor = FALSE)
plot(forelimb_lsr_pca) 
biplot(forelimb_lsr_pca) 
summary(forelimb_lsr_pca)
forelimb_lsr_pca_1_3 <- data.frame(forelimb_lsr_pca$scores[,1:3]) 


######### PhylogeneticEM ######### 
forelimb_lsr_pca_1_3_t <- t(forelimb_lsr_pca_1_3)

forelimb_lsr_pca_1_3_phyloEM <- PhyloEM(phylo = di2multi(tree_prune), 
                            Y_data = forelimb_lsr_pca_1_3_t, 
                            process = "scOU", 
                            random.root = TRUE, 
                            nbr_alpha = 20,             
                            parallel_alpha = TRUE,       
                            Ncores = 4)
save(forelimb_lsr_pca_1_3_phyloEM, file="forelimb_lsr_pca_1_3_phyloEM.Rdata")

plot_criterion(forelimb_lsr_pca_1_3_phyloEM)
plot(forelimb_lsr_pca_1_3_phyloEM, show.tip.label = TRUE, label_cex = .25)

all_params <- params_process(forelimb_lsr_pca_1_3_phyloEM, K = 9)
plot(forelimb_lsr_pca_1_3_phyloEM, params = all_params)

all_eq_sol <- equivalent_shifts(tree_prune, all_params)
plot(all_eq_sol, show_shifts_values = TRUE, shifts_cex = 0.5)



load("Rdata_hind/230531_data_hind.RData")
hindlimb_lsr <- as.matrix(select(data_hind, ends_with("_lsr_hind"))) #select log shape ratio traits

######### hindlimb Standard PCA ######### 
hindlimb_lsr_pca <- princomp(hindlimb_lsr, cor = FALSE)
plot(hindlimb_lsr_pca)
biplot(hindlimb_lsr_pca)
summary(hindlimb_lsr_pca)
hindlimb_lsr_pca_1_3 <- data.frame(hindlimb_lsr_pca$scores[,1:3])



######### PhylogeneticEM ######### 
hindlimb_lsr_pca_1_3_t <- t(hindlimb_lsr_pca_1_3)

hindlimb_lsr_pca_1_3_phyloEM <- PhyloEM(phylo = di2multi(tree_prune_hind), 
                                        Y_data = hindlimb_lsr_pca_1_3_t, 
                                        process = "scOU", 
                                        random.root = TRUE, 
                                        nbr_alpha = 20,             
                                        parallel_alpha = TRUE,       
                                        Ncores = 4)
save(hindlimb_lsr_pca_1_3_phyloEM, file="hindlimb_lsr_pca_1_3_phyloEM.Rdata")

plot_criterion(hindlimb_lsr_pca_1_3_phyloEM)
hindlimb_cols <- viridis(14)

plot(hindlimb_lsr_pca_1_3_phyloEM, show.tip.label = TRUE, label_cex = .25, automatic_colors=FALSE, color_characters= hindlimb_cols, color_edges = hindlimb_cols, color_shifts_regimes = TRUE)

all_params_hindlimb <- params_process(hindlimb_lsr_pca_1_3_phyloEM, K = 14)
plot(hindlimb_lsr_pca_1_3_phyloEM, params = all_params_hindlimb, show.tip.label = TRUE, label_cex = .25, gray_scale=TRUE)

all_eq_sol <- equivalent_shifts(tree_prune, all_params)
plot(all_eq_sol, show_shifts_values = TRUE, shifts_cex = 0.5)