#### Make character maps of locomotory regimes for hindlimb traits ####

#### libraries ####
library(tidyverse)
library(plyr)
library(plotrix)
library(geiger)
library(phytools)
library(phangorn)


#### read data ####
data_hind  <- read.csv("data_hind_loc3b.csv") #read species mean data
rownames(data_hind) <- data_hind[,1]

#### read multitree ####
multitree_hind <- read.nexus(file = "Upham_tree_batforelimbs_V2.nex")
multitree_hind <- lapply(multitree_hind, force.ultrametric)
multitree_hind <- lapply(multitree_hind, ladderize)

multitree_prune_hind <- list(length=length(multitree_hind))
for(i in 1:length(multitree_hind)){
  multitree_prune_hind[[i]] <- treedata(phy = multitree_hind[[i]], data = data_hind, warnings=TRUE)$phy
}
class(multitree_prune_hind) <- "multiPhylo"

#### obtain mcc tree from 1000 trees #### 
#tree_mcc <- maxCladeCred(multitree)
#write.tree(tree_mcc, file = "Upham_mcctree_batforelimbs_V2.tre")
tree_mcc_hind <- read.tree("Upham_mcctree_batforelimbs_V2.tre") #read mcc tree
tree_prune_hind <- treedata(tree_mcc_hind, data_hind, sort=TRUE)$phy #check that all species are in tree (nothing should appear)

data_hind <- data_hind[tree_prune_hind$tip.label,] #Sorts the rows of data so that they match the order of the tip labels. SUPER IMPORTANT


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

loc3_hind <- data_hind$Loc_mode_3groups  # changed column name to Loc_mode_3cgroups to avoid confusion with forelimb
names(loc3_hind) <- data_hind[,1]

loc3b_hind <- data_hind$Loc_mode_3bgroups
names(loc3b_hind) <- data_hind[,1]

loc3c_hind <- data_hind$Loc_mode_3cgroups
names(loc3c_hind) <- data_hind[,1]

loc4_hind <- data_hind$Loc_mode_4groups
names(loc4_hind) <- data_hind[,1]

#### make simmaps ####
multitree_prune_500_hind <- sample(multitree_prune_hind[1:1000],500) #subsample 500 of 1k  trees

### loc3c (generalists, arborealists, gliders+bats) ###
### all files are saved as just "loc3" bc we didn't make same mistake as forelimb
simmap_loc3_5k_hind <- make.simmap(multitree_prune_500_hind,loc3_hind, nsim=10) #10 simmap simulations across 500 trees
save(simmap_loc3_5k_hind, file="simmap_loc3_5k_hind.Rdata")
simmap_loc3_5k_s250_hind <- sample(simmap_loc3_5k_hind[1:5000],250) #subsample 250 of 5k simmap trees
save(simmap_loc3_5k_s250_hind, file="simmap_loc3_5k_s250_hind.Rdata")
simmap_loc3_mcc_s50_hind <- make.simmap(tree_prune_hind,loc3_hind, nsim=50) #10 simmap simulations across 1k trees
save(simmap_loc3_mcc_s50_hind, file="simmap_loc3_mcc_s50_hind.Rdata")

### loc3b (generalists, arb+gliders, bats) ###
simmap_loc3b_5k_hind <- make.simmap(multitree_prune_500_hind,loc3b_hind, nsim=10) #10 simmap simulations across 500 trees
save(simmap_loc3b_5k_hind, file="simmap_loc3b_5k_hind.Rdata")
simmap_loc3b_5k_s250_hind <- sample(simmap_loc3b_5k_hind[1:5000],250) #subsample 250 of 5k simmap trees
save(simmap_loc3b_5k_s250_hind, file="simmap_loc3b_5k_s250_hind.Rdata")

simmap_loc3b_mcc_s50_hind <- make.simmap(tree_prune_hind,loc3b_hind, nsim=50) #10 simmap simulations across 1k trees
save(simmap_loc3b_mcc_s50_hind, file="simmap_loc3b_mcc_s50_hind.Rdata")

### loc3c (generalists+arb, gliders, bats) ###
simmap_loc3c_5k_hind <- make.simmap(multitree_prune_500_hind,loc3c_hind, nsim=10) #10 simmap simulations across 500 trees
save(simmap_loc3c_5k_hind, file="simmap_loc3c_5k_hind.Rdata")
simmap_loc3c_5k_s250_hind <- sample(simmap_loc3c_5k_hind[1:5000],250) #subsample 250 of 5k simmap trees
save(simmap_loc3c_5k_s250_hind, file="simmap_loc3c_5k_s250_hind.Rdata")

simmap_loc3c_mcc_s50_hind <- make.simmap(tree_prune_hind,loc3c_hind, nsim=50) #10 simmap simulations across 1k trees
save(simmap_loc3c_mcc_s50_hind, file="simmap_loc3c_mcc_s50_hind.Rdata")

### loc 4  (generalists, arb, gliders, bats) ###
simmap_loc4_5k_hind <- make.simmap(multitree_prune_500_hind,loc4_hind, nsim=10) #10 simmap simulations across 500 trees
save(simmap_loc4_5k_hind, file="simmap_loc4_5k_hind.Rdata")
simmap_loc4_5k_s250_hind <- sample(simmap_loc4_5k_hind[1:5000],250) #subsample 250 of 5k simmap trees
save(simmap_loc4_5k_s250_hind, file="simmap_loc4_5k_s250_hind.Rdata")
simmap_loc4_mcc_s50_hind <- make.simmap(tree_prune_hind,loc4_hind, nsim=50) #10 simmap simulations across 1k trees
save(simmap_loc4_mcc_s50_hind, file="simmap_loc4_mcc_s50_hind.Rdata")