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

##### packages #####
library(mvMORPH)
library(tidyverse)
library(plyr)
library(plotrix)
library(geiger)
library(phytools)
library(phangorn)
options(scipen=999)

##### load data #####
load("Rdata_hind/230531_data_hind.RData")
load("Rdata_hind/simmap_hind/simmap_loc3_5k_s250_hind.Rdata")
load("Rdata_hind/simmap_hind/simmap_loc3b_5k_s250_hind.Rdata")
load("Rdata_hind/simmap_hind/simmap_loc3c_5k_s250_hind.Rdata")
load("Rdata_hind/simmap_hind/simmap_loc4_5k_s250_hind.Rdata")

##### variables #####
hindlimb_lsr <- as.matrix(select(data_hind, ends_with("_lsr_hind"))) #select log shape ratio traits

loc3_hind <- as.factor(data_hind$Loc_mode_3groups) #1_gen, 2_a, 3_g_fly
names(loc3_hind) <- data_hind[,1]

loc3b_hind <- as.factor(data_hind$Loc_mode_3bgroups) #1_gen, 2_a, 3_g_fly
names(loc3b_hind) <- data_hind[,1]

loc3c_hind <- as.factor(data_hind$Loc_mode_3cgroups) #1_gen, 2_a, 3_g_fly
names(loc3c_hind) <- data_hind[,1]

loc4_hind <- as.factor(data_hind$Loc_mode_4groups) #1_gen, 2_a, 3_g_fly
names(loc4_hind) <- data_hind[,1]


######### 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])


#### hindlimb mvMorph ####
# BM1 - traits follow brownian motion
hindlimb_BM1_pca <- list(length=length(simmap_loc4_5k_s250_hind))
for(i in 1:length(simmap_loc4_5k_s250_hind)){
  hindlimb_BM1_pca[[i]] <- mvBM(tree = simmap_loc4_5k_s250_hind[[i]], data = hindlimb_lsr_pca_1_3, model= "BM1", param= list(smean = TRUE))
  print(i)
}
save(hindlimb_BM1_pca, file="hindlimb_BM1_pca_1_3.Rdata")

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

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

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

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

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


### AIC model selection ####
load("Rdata_hind/mvMorph_models_hind/hindlimb_BM1_pca_1_3.Rdata")
load("Rdata_hind/mvMorph_models_hind/hindlimb_OU1_pca_1_3.Rdata")
load("Rdata_hind/mvMorph_models_hind/hindlimb_OUM_loc3_pca_1_3.Rdata")
load("Rdata_hind/mvMorph_models_hind/hindlimb_OUM_loc3b_pca_1_3.Rdata")
load("Rdata_hind/mvMorph_models_hind/hindlimb_OUM_loc3c_pca_1_3.Rdata")
load("Rdata_hind/mvMorph_models_hind/hindlimb_OUM_loc4_pca_1_3.Rdata")

hindlimb_PC_aic <- matrix(,length(hindlimb_OUM_loc4_pca),6,dimnames=list(c(1: length(hindlimb_OUM_loc4_pca)),c("mvBM1", "mvOU1", "mvOUM_loc3a","mvOUM_loc3b","mvOUM_loc3c","mvOUM_loc4")))
for(i in 1:length(hindlimb_OUM_loc4_pca)){
  hindlimb_PC_aic[i,1] <- hindlimb_BM1_pca[[i]]$AICc
  hindlimb_PC_aic[i,2] <- hindlimb_OU1_pca[[i]]$AICc
  hindlimb_PC_aic[i,3] <- hindlimb_OUM_loc3_pca[[i]]$AICc
  hindlimb_PC_aic[i,4] <- hindlimb_OUM_loc3b_pca[[i]]$AICc
  hindlimb_PC_aic[i,5] <- hindlimb_OUM_loc3c_pca[[i]]$AICc
  hindlimb_PC_aic[i,6] <- hindlimb_OUM_loc4_pca[[i]]$AICc
}
hindlimb_PC_aic

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