Source code

Full code for Figure 3 density plots, including Fig.3 E.

# load results

load("multifit_hum_fin.RData")
load("multifit_rad_fin.RData")
load("multifit_met_fin.RData")
load("multifit_phal_fin.RData")

## determinant 
sdet_hum <- sapply(multifit_hum, function(x) (det(x$sigma$Pinv)^(1/nrow(x$sigma$Pinv))))
sdet_rad <- sapply(multifit_rad, function(x) (det(x$sigma$Pinv)^(1/nrow(x$sigma$Pinv))))
sdet_met <- sapply(multifit_met, function(x) (det(x$sigma$Pinv)^(1/nrow(x$sigma$Pinv))))
sdet_phal <- sapply(multifit_phal, function(x) (det(x$sigma$Pinv)^(1/nrow(x$sigma$Pinv))))

det_traits <- cbind(sdet_hum,sdet_rad,sdet_met,sdet_phal)
colnames(det_traits) <- c("Humerus", "Radius", "Metacarpal", "Phalanx")
require(reshape2)
df_det<- melt(det_traits)

## trace
t_hum <- sapply(multifit_hum, function(x) sum(diag(x$sigma$Pinv)))
t_rad <- sapply(multifit_rad, function(x) sum(diag(x$sigma$Pinv)))
t_met <- sapply(multifit_met, function(x) sum(diag(x$sigma$Pinv)))
t_phal <- sapply(multifit_phal, function(x) sum(diag(x$sigma$Pinv)))

t_traits <- cbind(t_hum,t_rad,t_met,t_phal)
colnames(t_traits) <- c("Humerus", "Radius", "Metacarpal", "Phalanx")
df_t <- melt(t_traits)

## integration
library(evolqg)
e_hum <- sapply(multifit_hum, function(x) CalcEigenVar(cov2cor(x$sigma$Pinv), sd=TRUE))
e_rad <- sapply(multifit_rad, function(x) CalcEigenVar(cov2cor(x$sigma$Pinv), sd=TRUE))
e_met <- sapply(multifit_met, function(x) CalcEigenVar(cov2cor(x$sigma$Pinv), sd=TRUE))
e_phal <- sapply(multifit_phal, function(x) CalcEigenVar(cov2cor(x$sigma$Pinv), sd=TRUE))

i_traits <- cbind(e_hum,e_rad,e_met,e_phal)
colnames(i_traits) <- c("Humerus", "Radius", "Metacarpal", "Phalanx")
df_i <- melt(i_traits)

## rates
sigmas <- read.csv("sigmas.csv", sep=",", header=T)
sigmas <- sigmas[,-1]
med_sv <- apply(sigmas, 2, median)
med_sv_df <- as.data.frame(med_sv)
colnames(sigmas) <- c("Humerus","Radius","Metacarpal","Phalanx")
df_s <- melt(sigmas)
colnames(df_s) <- c("Bone", "Sig")

# original values:
load("results_raw.RData")
results

## 2) ridgeline to include in density plots

library(ggridges)
library(ggplot2)
library(viridis)
library(tidyverse)

cols<-setNames(c("#FFEB42", "#F7A24C", "#F163D1","#688CFD"),c("Humerus", "Radius", "Metacarpal", "Phalanx"))
library(pals)

resdf <- as.data.frame(results)

# Reorder following a precise order
p_det <- df_det %>%
 mutate(name = fct_relevel(Var2, 
                           "Phalanx", "Metacarpal", "Radius", "Humerus")) %>%
 ggplot(aes(x=value, y=name)) +
 theme_classic()+
 geom_density_ridges(alpha =0.6, aes(fill= Var2)) +
 scale_fill_manual(values =cols)+
 geom_vline(data=resdf, aes(xintercept=determinant, fill="black"), linetype = "dotted", size = 0.5)+
 xlab(NULL)+
 theme(legend.position="none")+
 ylab("determinant")

#include mean

p_t <- df_t %>%
 mutate(name = fct_relevel(Var2, 
                           "Phalanx", "Metacarpal", "Radius", "Humerus")) %>%
 ggplot(aes(x=value, y=name)) +
 theme_classic()+
 geom_density_ridges(alpha =0.6, aes(fill= Var2)) +
 scale_fill_manual(values =cols)+
 geom_vline(data=resdf, aes(xintercept=trace, fill="black"), linetype = "dotted", size = 0.5)+
 xlab(NULL)+
 theme(legend.position="none")+
 ylab("trace")

p_i <- df_i %>%
 mutate(name = fct_relevel(Var2, 
                           "Phalanx", "Metacarpal", "Radius", "Humerus")) %>%
 ggplot(aes(x=value, y=name)) +
 theme_classic()+
 geom_density_ridges(alpha =0.6, aes(fill= Var2)) +
 scale_fill_manual(values =cols)+
 geom_vline(data=resdf, aes(xintercept=integration, fill="black"), linetype = "dotted", size = 0.5)+
 xlab(NULL)+
 theme(legend.position="none")+
 ylab("integration")

p_s <- df_s %>%
 mutate(name = fct_relevel(Bone, 
                           "Phalanx", "Metacarpal", "Radius", "Humerus")) %>%
 ggplot(aes(x=Sig, y=Bone)) +
 theme_classic()+
 geom_density_ridges(alpha =0.6,aes(fill= Bone)) +
 scale_fill_manual(values =cols)+
 geom_vline(data=med_sv_df, aes(xintercept=med_sv, fill="black"), linetype = "dotted", size = 0.5)+
 xlab(NULL)+
 theme(legend.position="none")+
 ylab("stationary variance")

library(cowplot)
plot_grid(p_det, p_t, p_i, p_s, ncol=1, labels=c("B", "C", "D", "E"))