##### Get mean estimates for best extinct OUwie model for hdw (forelimb) trait #####
##### Saves mean optima results as csv #####

##### libraries #####
library(plotrix)

##### load data #####
load("Rdata_fore/230614_data_fore_ext.RData")
load("Rdata_fore/simmap_fore/simmap_loc3b_mcc_ext_s50.Rdata")
load("Rdata_fore/simmap_fore/simmap_loc3c_mcc_ext_s50.Rdata")
load("Rdata_fore/simmap_fore/simmap_loc3a_mcc_ext_s50.Rdata")
load("Rdata_fore/simmap_fore/simmap_loc4_mcc_ext_s50.Rdata")

##### load model results and variables #####
load("Rdata_fore/OUwie_models_fore_ext/hdw_lsr_OUM_loc3b_ext.Rdata")
options(scipen = 999)

species <- data$species
loc3b <- data$Loc_mode_3bgroups
hdw_lsr <- (as.numeric(as.character(data$hdw_lsr_fore)))
names(hdw_lsr) <- rownames(hdw_lsr)
hdw_lsr_loc3b <- data.frame(species, loc3b, hdw_lsr)

##### set up tables #####
# alpha constant
hdw_lsr_OUM_loc3b_ext_alpha <- matrix(,length(simmap_loc3b_mcc_ext_s50),1,dimnames=list(c(1: length(simmap_loc3b_mcc_ext_s50)),c("hdw_lsr_OUM_loc3b_ext_alpha")))
for(i in 1:length(simmap_loc3b_mcc_ext_s50)){
  hdw_lsr_OUM_loc3b_ext_alpha[i,1] <- hdw_lsr_OUM_loc3b_ext[[i]]$solution[1,1]
}
mean(hdw_lsr_OUM_loc3b_ext_alpha)
log(2)/mean(hdw_lsr_OUM_loc3b_ext_alpha)

# sigma constant
hdw_lsr_OUM_loc3b_ext_sigma <- matrix(,length(simmap_loc3b_mcc_ext_s50),1,dimnames=list(c(1: length(simmap_loc3b_mcc_ext_s50)),c("hdw_lsr_OUM_loc3b_ext_sigma")))
for(i in 1:length(simmap_loc3b_mcc_ext_s50)){
  hdw_lsr_OUM_loc3b_ext_sigma[i,1] <- hdw_lsr_OUM_loc3b_ext[[i]]$solution[2,1]
}
mean(hdw_lsr_OUM_loc3b_ext_sigma)

# optima regimes 
hdw_lsr_OUM_loc3b_ext_optima <- matrix(,length(simmap_loc3b_mcc_ext_s50),3, dimnames=list(c(1: length(simmap_loc3b_mcc_ext_s50)), c("1_gen", "2_a_g", "3_fly")))
for(i in 1:length(simmap_loc3b_mcc_ext_s50)){
  hdw_lsr_OUM_loc3b_ext_optima[i,1] <- hdw_lsr_OUM_loc3b_ext[[i]]$theta[1,1]
  hdw_lsr_OUM_loc3b_ext_optima[i,2] <- hdw_lsr_OUM_loc3b_ext[[i]]$theta[2,1]
  hdw_lsr_OUM_loc3b_ext_optima[i,3] <- hdw_lsr_OUM_loc3b_ext[[i]]$theta[3,1]
}
hdw_lsr_OUM_loc3b_ext_optima

mean_hdw_lsr_OUM_loc3b_ext_optima <- rbind(
  mean(hdw_lsr_OUM_loc3b_ext_optima[,1]),
  mean(hdw_lsr_OUM_loc3b_ext_optima[,2]),
  mean(hdw_lsr_OUM_loc3b_ext_optima[,3]))

CI_hdw_lsr_OUM_loc3b_ext_optima <- rbind(
  quantile(hdw_lsr_OUM_loc3b_ext_optima[,1], probs=c(0.025, 0.975)),
  quantile(hdw_lsr_OUM_loc3b_ext_optima[,2], probs=c(0.025, 0.975)),
  quantile(hdw_lsr_OUM_loc3b_ext_optima[,3], probs=c(0.025, 0.975)))

rownames(mean_hdw_lsr_OUM_loc3b_ext_optima) <- c("1_gen", "2_a_g", "3_fly")

##### save mean optima results as csv #####
mean_hdw_lsr_OUM_loc3b_ext_optima <- data.frame(mean_hdw_lsr_OUM_loc3b_ext_optima,CI_hdw_lsr_OUM_loc3b_ext_optima)
colnames(mean_hdw_lsr_OUM_loc3b_ext_optima) <- c("OUM_optima", "L95", "U95")
mean_hdw_lsr_OUM_loc3b_ext_optima <- mean_hdw_lsr_OUM_loc3b_ext_optima-10
mean_hdw_lsr_OUM_loc3b_ext_optima
write.csv(mean_hdw_lsr_OUM_loc3b_ext_optima, file="output_fore/OUwie_mean_optima_ext/mean_hdw_lsr_OUM_loc3b_ext_optima_ext.csv")