##### Bootstrapping for hdw (forelimb) trait ##### ##### This is bootstrapping for loc3b since that was the best fit ##### ##### model, but you can switch out the variables for loc4, etc. ##### ##### Outputs optima means and CIs as csv and saves plots as pdfs ##### ##### libraries ##### library(OUwie) library(plotrix) ##### load data ##### load("230531_data_fore.RData") load("simmap_loc3b_5k_s250.Rdata") load("hdw_lsr_OUM_loc3b.Rdata") ##### variables ##### species <- data$?..species loc3b <- data$Loc_mode_3bgroups hdw_lsr <- (as.numeric(as.character(data$hdw_lsr_fore))) names(hdw_lsr) <- rownames(hdw_lsr) #### hindlimb dataframes for OUwie #### hdw_lsr_loc3b <- data.frame(species, loc3b, hdw_lsr) options(scipen = 999) # alpha constant hdw_lsr_OUM_loc3b_alpha <- matrix(,length(simmap_loc3b_5k_s250),1,dimnames=list(c(1: length(simmap_loc3b_5k_s250)),c("hdw_lsr_OUM_loc3b_alpha"))) for(i in 1:length(simmap_loc3b_5k_s250)){ hdw_lsr_OUM_loc3b_alpha[i,1] <- hdw_lsr_OUM_loc3b[[i]]$solution[1,1] } mean(hdw_lsr_OUM_loc3b_alpha) log(2)/mean(hdw_lsr_OUM_loc3b_alpha) # sigma constant hdw_lsr_OUM_loc3b_sigma <- matrix(,length(simmap_loc3b_5k_s250),1,dimnames=list(c(1: length(simmap_loc3b_5k_s250)),c("hdw_lsr_OUM_loc3b_sigma"))) for(i in 1:length(simmap_loc3b_5k_s250)){ hdw_lsr_OUM_loc3b_sigma[i,1] <- hdw_lsr_OUM_loc3b[[i]]$solution[2,1] } mean(hdw_lsr_OUM_loc3b_sigma) # optima regimes hdw_lsr_OUM_loc3b_optima <- matrix(,length(simmap_loc3b_5k_s250),3, dimnames=list(c(1: length(simmap_loc3b_5k_s250)), c("1_gen", "2_a_g", "3_fly"))) for(i in 1:length(simmap_loc3b_5k_s250)){ hdw_lsr_OUM_loc3b_optima[i,1] <- hdw_lsr_OUM_loc3b[[i]]$theta[1,1] hdw_lsr_OUM_loc3b_optima[i,2] <- hdw_lsr_OUM_loc3b[[i]]$theta[2,1] hdw_lsr_OUM_loc3b_optima[i,3] <- hdw_lsr_OUM_loc3b[[i]]$theta[3,1] } hdw_lsr_OUM_loc3b_optima mean_hdw_lsr_OUM_loc3b_optima <- rbind( mean(hdw_lsr_OUM_loc3b_optima[,1]), mean(hdw_lsr_OUM_loc3b_optima[,2]), mean(hdw_lsr_OUM_loc3b_optima[,3])) CI_hdw_lsr_OUM_loc3b_optima <- rbind( quantile(hdw_lsr_OUM_loc3b_optima[,1], probs=c(0.025, 0.975)), quantile(hdw_lsr_OUM_loc3b_optima[,2], probs=c(0.025, 0.975)), quantile(hdw_lsr_OUM_loc3b_optima[,3], probs=c(0.025, 0.975))) rownames(mean_hdw_lsr_OUM_loc3b_optima) <- c("1_gen", "2_a_g", "3_fly") mean_hdw_lsr_OUM_loc3b_optima <- data.frame(mean_hdw_lsr_OUM_loc3b_optima,CI_hdw_lsr_OUM_loc3b_optima) colnames(mean_hdw_lsr_OUM_loc3b_optima) <- c("OUM_optima", "L95", "U95") mean_hdw_lsr_OUM_loc3b_optima <- mean_hdw_lsr_OUM_loc3b_optima-10 mean_hdw_lsr_OUM_loc3b_optima #### Bootstrapping for hdw_lsr_OUM_loc3b model ##### simmap_loc3b_5k_s250_10 <- sample(simmap_loc3b_5k_s250[1:250],10) hdw_lsr_OUM_loc3b_boot <- list(length=length(simmap_loc3b_5k_s250_10)) for(i in 1:length(simmap_loc3b_5k_s250_10)){ hdw_lsr_OUM_loc3b_boot[[i]] <- OUwie.boot(simmap_loc3b_5k_s250_10[[i]], hdw_lsr_loc3b, model="OUM",simmap.tree=TRUE, root.station=TRUE, nboot=100, mserr="none", alpha=c(mean(hdw_lsr_OUM_loc3b_alpha), mean(hdw_lsr_OUM_loc3b_alpha), mean(hdw_lsr_OUM_loc3b_alpha)), sigma.sq=c(mean(hdw_lsr_OUM_loc3b_sigma), mean(hdw_lsr_OUM_loc3b_sigma), mean(hdw_lsr_OUM_loc3b_sigma)), theta=c(mean(hdw_lsr_OUM_loc3b_optima[,1]), mean(hdw_lsr_OUM_loc3b_optima[,2]), mean(hdw_lsr_OUM_loc3b_optima[,3])), theta0=mean(hdw_lsr_OUM_loc3b_optima[,1])) print(i) } save(hdw_lsr_OUM_loc3b_boot, file="hdw_lsr_OUM_loc3b_boot.Rdata") ##### results ##### load("Rdata_fore/OUwie_boot_fore/hdw_lsr_OUM_loc3b_boot.Rdata") ### hdw_lsr_OUM_loc3b optima ### hdw_lsr_OUM_loc3b_boot_opt_2_gen <- list() for(i in 1:length(hdw_lsr_OUM_loc3b_boot)){ hdw_lsr_OUM_loc3b_boot_opt_2_gen[[i]] <- c(hdw_lsr_OUM_loc3b_boot[[i]][,7])} hdw_lsr_OUM_loc3b_boot_opt_4_a_g <- list() for(i in 1:length(hdw_lsr_OUM_loc3b_boot)){ hdw_lsr_OUM_loc3b_boot_opt_4_a_g[[i]] <- c(hdw_lsr_OUM_loc3b_boot[[i]][,8])} hdw_lsr_OUM_loc3b_boot_opt_5_fly <- list() for(i in 1:length(hdw_lsr_OUM_loc3b_boot)){ hdw_lsr_OUM_loc3b_boot_opt_5_fly[[i]] <- c(hdw_lsr_OUM_loc3b_boot[[i]][,9])} #confidence interval df mean_hdw_lsr_OUM_loc3b_boot_opt <- rbind( mean(unlist(hdw_lsr_OUM_loc3b_boot_opt_2_gen)), mean(unlist(hdw_lsr_OUM_loc3b_boot_opt_4_a_g)), mean(unlist(hdw_lsr_OUM_loc3b_boot_opt_5_fly))) mean_hdw_lsr_OUM_loc3b_boot_opt <- mean_hdw_lsr_OUM_loc3b_boot_opt-10 CI_hdw_lsr_OUM_loc3b_boot_opt <- rbind( quantile(unlist(hdw_lsr_OUM_loc3b_boot_opt_2_gen), probs=c(0.025, 0.975)), quantile(unlist(hdw_lsr_OUM_loc3b_boot_opt_4_a_g), probs=c(0.025, 0.975)), quantile(unlist(hdw_lsr_OUM_loc3b_boot_opt_5_fly), probs=c(0.025, 0.975))) CI_hdw_lsr_OUM_loc3b_boot_opt <- CI_hdw_lsr_OUM_loc3b_boot_opt-10 hdw_lsr_OUM_loc3b_boot_opt <- data.frame(mean_hdw_lsr_OUM_loc3b_optima$OUM_optima,CI_hdw_lsr_OUM_loc3b_boot_opt) colnames(hdw_lsr_OUM_loc3b_boot_opt) <- c("OUM_optima_mean", "L95", "U95") rownames(hdw_lsr_OUM_loc3b_boot_opt) <- c("2_gen", "4_a_g", "5_fly") hdw_lsr_OUM_loc3b_boot_opt write.csv(hdw_lsr_OUM_loc3b_boot_opt, file="hdw_lsr_OUM_loc3b_boot_opt.csv") ##### plot optima and confidence intervals ##### pdf(file="hdw.pdf", height=5, width=4) par(mar=c(1.5,5,1.5,5)+1) cols_loc3b <- c("#E6C060", "#337267","#AB97AF") plot(mean_hdw_lsr_OUM_loc3b_optima$OUM_optima, ylim = c(min(mean_hdw_lsr_OUM_loc3b_optima$OUM_optima-(hdw_lsr_OUM_loc3b_boot_opt$OUM_optima_mean-hdw_lsr_OUM_loc3b_boot_opt$L95)), max(hdw_lsr_OUM_loc3b_boot_opt$U95-hdw_lsr_OUM_loc3b_boot_opt$OUM_optima_mean+mean_hdw_lsr_OUM_loc3b_optima$OUM_optima)), xlim=c(0.5,3.5), cex=2.7, pch=20, col=cols_loc3b, xlab=" ", ylab=" ", xaxt = "n", las=1, cex.axis=1.3) title(ylab="humerus distal width optima", line=3.5, cex.lab=1.6) axis(1, xaxp=c(1, 3, 2), labels=F) plotCI(mean_hdw_lsr_OUM_loc3b_optima$OUM_optima,y=NULL, uiw=hdw_lsr_OUM_loc3b_boot_opt$U95-hdw_lsr_OUM_loc3b_boot_opt$OUM_optima_mean, liw=hdw_lsr_OUM_loc3b_boot_opt$OUM_optima_mean-hdw_lsr_OUM_loc3b_boot_opt$L95, err="y", pch=20, slty=1, col=cols_loc3b, scol = cols_loc3b, add=TRUE, lwd=7) dev.off()