##### OUwie models for hdw (forelimb) trait #####
##### Runs loc4, loc3 (same as loc3a), loc3b, and loc3c models #####
##### Saves models Rdata and outputs csv with AICw scores #####

##### libraries #####
library(OUwie)
options(scipen=999)

##### load data #####
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")

##### variables #####
species <- data$?..species
loc3 <- data$Loc_mode_3groups
loc3b <- data$Loc_mode_3bgroups
loc3c <- data$Loc_mode_3cgroups
loc4 <- data$Loc_mode_4groups

### hdw ###
hdw_lsr <- (as.numeric(as.character(data$hdw_lsr)))
names(hdw_lsr) <- rownames(hdw_lsr)

### forelimb dataframes for OUwie ###
hdw_lsr_loc3 <- data.frame(species, loc3, (hdw_lsr+10))
hdw_lsr_loc3b <- data.frame(species, loc3b, (hdw_lsr+10))
hdw_lsr_loc3c <- data.frame(species, loc3c, (hdw_lsr+10))
hdw_lsr_loc4 <- data.frame(species, loc4, (hdw_lsr+10))


##### OUwie models for hdw_lsr #####
### BM model ###
hdw_lsr_BM <- list(length=length(simmap_loc4_5k_s250))
for(i in 1:length(simmap_loc4_5k_s250)){
  hdw_lsr_BM[[i]] <- OUwie(simmap_loc4_5k_s250[[i]], hdw_lsr_loc4, model="BM1", simmap.tree=TRUE,root.station=FALSE, mserr = "none", algorithm = "three.point", diagn=T)
  print(i)
}
save(hdw_lsr_BM, file="hdw_lsr_models_v6/hdw_lsr_BM.Rdata")


### OU model ###
hdw_lsr_OU <- list(length=length(simmap_loc4_5k_s250))
for(i in 1:length(simmap_loc4_5k_s250)){
  hdw_lsr_OU[[i]] <- OUwie(simmap_loc4_5k_s250[[i]], hdw_lsr_loc4, model="OU1", simmap.tree=TRUE,root.station=FALSE, mserr = "none", algorithm = "three.point", diagn=T )
  print(i)
}
save(hdw_lsr_OU, file="hdw_lsr_models_v6/hdw_lsr_OU.Rdata")


#### loc3 OUM Models ####
hdw_lsr_OUM_loc3 <- list(length=length(simmap_loc3_5k_s250))
for(i in 1:length(simmap_loc3_5k_s250)){
  hdw_lsr_OUM_loc3[[i]] <- OUwie(simmap_loc3_5k_s250[[i]], hdw_lsr_loc3, model="OUM", simmap.tree=TRUE,root.station=FALSE, mserr = "none", algorithm = "three.point", diagn=T)
  print(i)
}
save(hdw_lsr_OUM_loc3, file="hdw_lsr_OUM_loc3.Rdata")


#### loc3b OUM Models ####
hdw_lsr_OUM_loc3b <- list(length=length(simmap_loc3b_5k_s250))
for(i in 1:length(simmap_loc3b_5k_s250)){
  hdw_lsr_OUM_loc3b[[i]] <- OUwie(simmap_loc3b_5k_s250[[i]], hdw_lsr_loc3b, model="OUM", simmap.tree=TRUE,root.station=FALSE, mserr = "none", algorithm = "three.point", diagn=T)
  print(i)
}
save(hdw_lsr_OUM_loc3b, file="hdw_lsr_OUM_loc3b.Rdata")


#### loc3c OUM Models ####
hdw_lsr_OUM_loc3c <- list(length=length(simmap_loc3c_5k_s250))
for(i in 1:length(simmap_loc3c_5k_s250)){
  hdw_lsr_OUM_loc3c[[i]] <- OUwie(simmap_loc3c_5k_s250[[i]], hdw_lsr_loc3c, model="OUM", simmap.tree=TRUE,root.station=FALSE, mserr = "none", algorithm = "three.point", diagn=T)
  print(i)
}
save(hdw_lsr_OUM_loc3c, file="hdw_lsr_OUM_loc3c.Rdata")


#### loc4 OUM Models ####
hdw_lsr_OUM_loc4 <- list(length=length(simmap_loc4_5k_s250))
for(i in 1:length(simmap_loc4_5k_s250)){
  hdw_lsr_OUM_loc4[[i]] <- OUwie(simmap_loc4_5k_s250[[i]], hdw_lsr_loc4, model="OUM", simmap.tree=TRUE,root.station=FALSE, mserr = "none", algorithm = "three.point", diagn=T)
  print(i)
}
save(hdw_lsr_OUM_loc4, file="hdw_lsr_OUM_loc4.Rdata")


#### load all OUwie models (if needed) ####
load("Rdata_fore/OUwie_models_fore/hdw_lsr_OUM_loc3.Rdata")
load("Rdata_fore/OUwie_models_fore/hdw_lsr_OUM_loc3b.Rdata")
load("Rdata_fore/OUwie_models_fore/hdw_lsr_OUM_loc3c.Rdata")
load("Rdata_fore/OUwie_models_fore/hdw_lsr_OUM_loc4.Rdata")
load("Rdata_fore/OUwie_models_fore/hdw_lsr_OU.Rdata")
load("Rdata_fore/OUwie_models_fore/hdw_lsr_BM.Rdata")

hdw_lsr_aic.table <- matrix(,length(hdw_lsr_BM),6,dimnames=list(c(1: length(hdw_lsr_BM)),c("BM1", "OU1", "OUM_loc3", "OUM_loc3b", "OUM_loc3c", "OUM_loc4" )))
for(i in 1:length(hdw_lsr_BM)){
  hdw_lsr_aic.table[i,1] <- hdw_lsr_BM[[i]]$AICc
  hdw_lsr_aic.table[i,2] <- hdw_lsr_OU[[i]]$AICc
  hdw_lsr_aic.table[i,3] <- hdw_lsr_OUM_loc3[[i]]$AICc
  hdw_lsr_aic.table[i,4] <- hdw_lsr_OUM_loc3b[[i]]$AICc
  hdw_lsr_aic.table[i,5] <- hdw_lsr_OUM_loc3c[[i]]$AICc
  hdw_lsr_aic.table[i,6] <- hdw_lsr_OUM_loc4[[i]]$AICc
}
hdw_lsr_aic.table

hdw_lsr_aic.table_aicw <- 
  geiger::aicw(c(
    mean(hdw_lsr_aic.table[,1]),
    mean(hdw_lsr_aic.table[,2]),
    mean(hdw_lsr_aic.table[,3]),
    mean(hdw_lsr_aic.table[,4]),
    mean(hdw_lsr_aic.table[,5]),
    mean(hdw_lsr_aic.table[,6])))
rownames(hdw_lsr_aic.table_aicw) <- c("BM1", "OU1",  "OUM_loc3", "OUM_loc3b", "OUM_loc3c", "OUM_loc4")
hdw_lsr_aic.table_aicw
write.csv(hdw_lsr_aic.table_aicw, file="hdw_lsr_aicw_table.csv")