##### MvMorph multivariate models with PC1-3 only #####
##### extant-only dataset #####

##### packages #####
library(tidyverse)
library(plyr)
library(plotrix)
library(geiger)
library(phytools)
library(mvMORPH)
library(tidyselect)

options(scipen=999)

load("Rdata_fore/230531_data_fore.RData")
load("Rdata_fore/simmap_fore/simmap_loc3_5k_s250.Rdata")
load("Rdata_fore/simmap_fore/simmap_loc3b_5k_s250.Rdata")
load("Rdata_fore/simmap_fore/simmap_loc3c_5k_s250.Rdata")
load("Rdata_fore/simmap_fore/simmap_loc4_5k_s250.Rdata")


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

loc4 <- as.factor(data$Loc_mode_4groups) #1_gen, 2_a, 3_g, 4_fly
names(loc4) <- data[,1]

loc3 <- as.factor(data$Loc_mode_3groups) #1_gen, 2_a, 3_g, 4_fly
names(loc3) <- data[,1]

loc3b <- as.factor(data$Loc_mode_3bgroups) #1_gen, 2_a, 3_g, 4_fly
names(loc3b) <- data[,1]

loc3c <- as.factor(data$Loc_mode_3cgroups) #1_gen, 2_a, 3_g, 4_fly
names(loc3c) <- data[,1]

######### 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]) 
# check loadings: forelimb_lsr_pca$loadings
save(forelimb_lsr_pca_1_3, file="forelimb_lsr_pca_1_3.Rdata")

#### plot PCA for visuzalization #####
cols_loc4 <- setNames(c("#E6C060", #generalist (white)
                                 "#419D78", #arboreal
                                 "#AC97AF", #gliders (tan)
                                 "#274754"), #flying (navy)
                                 levels(loc4))
points_loc4 = c(22, 21, 23, 24)

# forelimb PCA
plot(forelimb_lsr_pca_1_3[,1], forelimb_lsr_pca_1_3[,2], bg = cols_loc4[loc4], pch = points_loc4[loc4], cex = 1.5, las =1)
for(jj in 1:length(levels(loc4))){
  ii=levels(loc4)[jj]
  kk=chull(forelimb_lsr_pca_1_3[loc4==ii,1:2])
  lines(forelimb_lsr_pca_1_3[loc4==ii,1][c(kk, kk[1])],
        forelimb_lsr_pca_1_3[loc4==ii,2][c(kk, kk[1])], col="black")
}
text(forelimb_lsr_pca_1_5[,1], forelimb_lsr_pca_1_5[,2], labels = rownames(data), pos = 4, cex = 0.5)

#### forelimb mvMorph ####
# BM1 - traits follow brownian motion
forelimb_BM1_pca <- list(length=length(simmap_loc4_5k_s250))
for(i in 1:length(simmap_loc4_5k_s250)){
  forelimb_BM1_pca[[i]] <- mvBM(tree = simmap_loc4_5k_s250[[i]], data = forelimb_lsr_pca_1_3, model= "BM1", param= list(smean = TRUE))
  print(i)
}
save(forelimb_BM1_pca, file="forelimb_BM1_pca_1_3.Rdata")

# OU1
forelimb_OU1_pca <- list(length=length(simmap_loc4_5k_s250))
for(i in 1:length(simmap_loc4_5k_s250)){
  forelimb_OU1_pca[[i]] <- mvOU(tree = simmap_loc4_5k_s250[[i]], data = forelimb_lsr_pca_1_3, model= "OU1", param= list(root = "stationary"))
  print(i)
}
save(forelimb_OU1_pca, file="forelimb_OU1_pca_1_3.Rdata")

# OUM_loc3
forelimb_OUM_loc3_pca <- list(length=length(simmap_loc3_5k_s250))
for(i in 1:length(simmap_loc3_5k_s250)){
  forelimb_OUM_loc3_pca[[i]] <- mvOU(tree = simmap_loc3_5k_s250[[i]], data = forelimb_lsr_pca_1_3, model= "OUM", param= list(root = "stationary"))
  print(i)
}
save(forelimb_OUM_loc3_pca, file="forelimb_OUM_loc3_pca_1_3.Rdata")

# OUM_loc3b
forelimb_OUM_loc3b_pca <- list(length=length(simmap_loc3b_5k_s250))
for(i in 1:length(simmap_loc3b_5k_s250)){
  forelimb_OUM_loc3b_pca[[i]] <- mvOU(tree = simmap_loc3b_5k_s250[[i]], data = forelimb_lsr_pca_1_3, model= "OUM", param= list(root = "stationary"))
  print(i)
}
save(forelimb_OUM_loc3b_pca, file="forelimb_OUM_loc3b_pca_1_3.Rdata")

# OUM_loc3c
forelimb_OUM_loc3c_pca <- list(length=length(simmap_loc3c_5k_s250))
for(i in 1:length(simmap_loc3c_5k_s250)){
  forelimb_OUM_loc3c_pca[[i]] <- mvOU(tree = simmap_loc3c_5k_s250[[i]], data = forelimb_lsr_pca_1_3, model= "OUM", param= list(root = "stationary"))
  print(i)
}
save(forelimb_OUM_loc3c_pca, file="forelimb_OUM_loc3c_pca_1_3.Rdata")

# OUM_loc4
forelimb_OUM_loc4_pca <- list(length=length(simmap_loc4_5k_s250))
for(i in 1:length(simmap_loc4_5k_s250)){
  forelimb_OUM_loc4_pca[[i]] <- mvOU(tree = simmap_loc4_5k_s250[[i]], data = forelimb_lsr_pca_1_3, model= "OUM", param= list(root = "stationary"))
  print(i)
}
save(forelimb_OUM_loc4_pca, file="forelimb_OUM_loc4_pca_1_3.Rdata")

##### AIC table for forelimb models #####
load("Rdata_fore/mvMorph_models_fore/forelimb_BM1_pca_1_3.Rdata")
load("Rdata_fore/mvMorph_models_fore/forelimb_OU1_pca_1_3.Rdata")
load("Rdata_fore/mvMorph_models_fore/forelimb_OUM_loc3_pca_1_3.RData")
load("Rdata_fore/mvMorph_models_fore/forelimb_OUM_loc3b_pca_1_3.RData")
load("Rdata_fore/mvMorph_models_fore/forelimb_OUM_loc3c_pca_1_3.RData")
load("Rdata_fore/mvMorph_models_fore/forelimb_OUM_loc4_pca_1_3.RData")

forelimb_PC_aic.table <- matrix(,length(forelimb_OUM_loc4_pca),6,dimnames=list(c(1: length(forelimb_OUM_loc4_pca)),c("mvBM1", "mvOU1","mvOUM_loc3","mvOUM_loc3b","mvOUM_loc3c", "mvOUM_loc4")))
for(i in 1:length(forelimb_OUM_loc4_pca)){
  forelimb_PC_aic.table[i,1] <- forelimb_BM1_pca[[i]]$AICc
  forelimb_PC_aic.table[i,2] <- forelimb_OU1_pca[[i]]$AICc
  forelimb_PC_aic.table[i,3] <- forelimb_OUM_loc3_pca[[i]]$AICc
  forelimb_PC_aic.table[i,4] <- forelimb_OUM_loc3b_pca[[i]]$AICc
  forelimb_PC_aic.table[i,5] <- forelimb_OUM_loc3c_pca[[i]]$AICc
  forelimb_PC_aic.table[i,6] <- forelimb_OUM_loc4_pca[[i]]$AICc
}
forelimb_PC_aic.table

forelimb_PC_aicw <- 
  geiger::aicw(c(
    mean(forelimb_PC_aic.table[,1]),
    mean(forelimb_PC_aic.table[,2]),
    mean(forelimb_PC_aic.table[,3]),
    mean(forelimb_PC_aic.table[,4]),
    mean(forelimb_PC_aic.table[,5]),
    mean(forelimb_PC_aic.table[,6])
  ))
rownames(forelimb_PC_aicw) <- c("mvBM1", "mvOU1", "mvOUM_loc3","mvOUM_loc3b","mvOUM_loc3c","mvOUM_loc4")
forelimb_PC_aicw
write.csv(forelimb_PC_aicw, file="forelimb_PC1-3_aicw.csv")