---
title: "Problems with using comparative analyses of avian brain size to test hypotheses of cognitive evolution"
script author: RH
date created: 13th may 2021
output:
  html_document:
    toc: yes
    toc_depth: 4
  'html_document:': default
---
```{r include = FALSE}
# load libs
library("ape")
library("caper")
library("phytools")
library("corHMM")
library("plyr")
library("car")
library("dplyr")
library("ggplot2")
library("tidyverse")
library("RColorBrewer")
library("scales")
library("ggthemes")
library("gridExtra")

# load data
brain_size<- read.csv("brain_size.csv", header=TRUE)
body_size<- read.csv("body_size.csv", header=TRUE)
corvides<- read.csv("corvides.csv", header=TRUE)
tree<-read.nexus("corvides.nex")
tree2<-multi2di(tree)
```

# variation in brain and body size
## brain size
```{r include = FALSE}
nrow(brain_size)

# subset data to 2 or more datapoints per species (sex-averaged)
brain_size_ave <- brain_size
brain_size_ave$count <- rowSums( !is.na( brain_size_ave[,c(4:9)]))
brain_size_ave2 <- brain_size_ave %>%
  dplyr::filter(count > 1)
nrow(brain_size_ave2)

# remove sex-specific measure
brain_size_ave2 <- brain_size_ave2 %>%
  dplyr::select(! c("GarciaPenaMLog","GarciaPenaFLog"))
View(brain_size_ave2)
brain_size_ave3 <- gather(brain_size_ave2, "dataset", "brain_volume", c(2,3,4,5,6,7))
brain_size_ave3 <- brain_size_ave3 %>%
  arrange(brain_volume)
brain_size_ave3$Species <- as.character(brain_size_ave3$Species)
brain_size_ave3$Species <- factor((brain_size_ave3$Species), levels=unique(brain_size_ave3$Species))
brain_size_ave3$dataset <- as.factor(brain_size_ave3$dataset)

# remove rows w NA brain volume
brain_size_ave4 <- brain_size_ave3 %>%
  filter(!is.na(brain_volume))
summary(brain_size_ave4$dataset)

for(i in 1:nrow(brain_size_ave4)){
  if(brain_size_ave4$dataset[i] == "SayolLog"){
    brain_size_ave4$`dataset number`[i] <- "dataset 1"}
  if(brain_size_ave4$dataset[i] == "IwaniukNelsonLog"){
    brain_size_ave4$`dataset number`[i] <- "dataset 2"}
  if(brain_size_ave4$dataset[i] == "GarciaPenaAveLog"){
    brain_size_ave4$`dataset number`[i] <- "dataset 3"}
  if(brain_size_ave4$dataset[i] == "IwaniukArnoldLog"){
    brain_size_ave4$`dataset number`[i] <- "dataset 4"}
  if(brain_size_ave4$dataset[i] == "IwaniukDeanNelsonLog"){
    brain_size_ave4$`dataset number`[i] <- "dataset 5"}
  if(brain_size_ave4$dataset[i] == "CorfieldLog"){
    brain_size_ave4$`dataset number`[i] <- "dataset 6"}
}

brain_size_plot <- ggplot(brain_size_ave4, aes(x = brain_volume,Species))+
  scale_colour_brewer(palette = "Paired")+
  geom_line(aes(),size = .1, alpha =.3)+
  geom_point(aes(colour = `dataset number`, alpha = .2))+
  labs(tag="a")+
  xlab("log brain volume (ml)")+
  ylab("species")+
  expand_limits(y= c(-100, length(levels(brain_size_ave4$Species)) + 100))+
  theme_classic()+
  theme(
   axis.text.y = element_blank(),
   axis.ticks = element_blank(),
   legend.position = "none")
brain_size_plot
```

###  quantifying variation
```{r echo = FALSE}
brain_size_ave2$min <- (apply(brain_size_ave2[,2:7], 1, min, na.rm = TRUE))
brain_size_ave2$max <- (apply(brain_size_ave2[,2:7], 1, max, na.rm = TRUE))
brain_size_ave2$diff <- (brain_size_ave2$max) - (brain_size_ave2$min)
View(brain_size_ave2) # which have largest diff? Pyrrhura_frontalis, Acanthorhynchus_tenuirostris, Tinamus_major

brain_size_ave4 <- brain_size_ave4 %>%
  arrange(brain_volume)
brain_size_ave_all <- gather(brain_size_ave, "dataset", "brain_volume", c(4:9))
brain_size_ave_all <- brain_size_ave_all %>%
  dplyr::select(Species, dataset, brain_volume)

# remove rows w NA brain volume
brain_size_ave_all1 <- brain_size_ave_all %>%
  filter(!is.na(brain_volume))
brain_size_ave_all2 <- brain_size_ave_all1 %>%
  arrange(brain_volume)

#Pyrrhura_frontalis
which(brain_size_ave_all2$Species %in% "Pyrrhura_frontalis")
pf <- brain_size_ave_all2[c(1574:1862), ]
pf <- pf %>%
  dplyr::distinct(Species)
nrow(pf) # 159 species
(nrow(pf)/nrow(brain_size_ave))*100 # 10.79%

#Acanthorhynchus_tenuirostris
which(brain_size_ave_all2$Species %in% "Acanthorhynchus_tenuirostris")
at <- brain_size_ave_all2[c(213:473), ]
at <- at %>%
  dplyr::distinct(Species) # 150 distinct species
nrow(at) 
(nrow(at)/nrow(brain_size_ave))*100 # overlaps with 10.18%

#Tinamus_major
which(brain_size_ave_all2$Species %in% "Tinamus_major")
tm <- brain_size_ave_all2[c(1627:1885), ]
tm <- tm %>%
  dplyr::distinct(Species)
nrow(tm) #141
(nrow(tm)/nrow(brain_size_ave))*100 # 9.57%
```

## body size
```{r include = FALSE}
nrow(body_size)

# subset data to 2 or more datapoints per species (sex-averaged)
body_size$count <- rowSums( !is.na( body_size[,2:12]))

body_size2 <- body_size %>%
  dplyr::filter(count > 1)
nrow(body_size2) 
body_size3 <- gather(body_size2, "dataset", "body_mass", 2:12)

# remove rows NA body mass
body_size4 <- body_size3 %>%
  filter(!is.na(body_mass))

for(i in 1:nrow(body_size4)){
  if(body_size4$dataset[i] == "CorfieldLog"){
    body_size4$`dataset number`[i] <- "dataset 1"}
  if(body_size4$dataset[i] == "FristoeLog"){
    body_size4$`dataset number`[i] <- "dataset 2"}
  if(body_size4$dataset[i] == "GaramszegiLog"){
    body_size4$`dataset number`[i] <- "dataset 3"}
  if(body_size4$dataset[i] == "Iwaniuk2004Log"){
    body_size4$`dataset number`[i] <- "dataset 4"}
  if(body_size4$dataset[i] == "IwaniukDeanNelsonLog"){
    body_size4$`dataset number`[i] <- "dataset 5"}
    if(body_size4$dataset[i] == "IwaniukArnoldLog"){
    body_size4$`dataset number`[i] <- "dataset 6"}
  if(body_size4$dataset[i] == "LendvaiLog"){
    body_size4$`dataset number`[i] <- "dataset 7"}
  if(body_size4$dataset[i] == "Minias1Log"){
    body_size4$`dataset number`[i] <- "dataset 8"}
  if(body_size4$dataset[i] == "Minias2Log"){
    body_size4$`dataset number`[i] <- "dataset 9"}
  if(body_size4$dataset[i] == "SayolLog"){
    body_size4$`dataset number`[i] <- "dataset 10"}
  if(body_size4$dataset[i] == "SolLog"){
    body_size4$`dataset number`[i] <- "dataset 11"}
}

body_size4$`dataset number` <- as.factor(body_size4$`dataset number`)
body_size4$body_mass <- as.numeric((body_size4$body_mass))
body_size4 <- body_size4 %>%
  arrange(body_mass)
body_size4$Species <- as.character(body_size4$Species)
body_size4$Species <- factor((body_size4$Species), levels=unique(body_size4$Species))

body_size_plot <- ggplot(body_size4, aes(x = body_mass,y = Species))+
  scale_colour_brewer(palette = "Paired")+
  geom_line(aes(),size = .1, alpha =.3)+
  geom_point(aes(colour = `dataset number`, alpha = .2))+
  labs(tag="b")+
  xlab("log body mass (g)")+
  ylab("species")+
  expand_limits(y= c(-100, length(levels(body_size4$Species))+ 100))+
  theme_classic()+
  theme(
   axis.text.y = element_blank(),
   axis.ticks = element_blank(),
   legend.position = "none")
body_size_plot
```

###  quantifying variation
```{r echo = FALSE}
body_size2$min <- as.numeric((apply(body_size2[,2:12], 1, min, na.rm = TRUE)))
body_size2$max <- as.numeric((apply(body_size2[,2:12], 1, max, na.rm = TRUE)))
body_size2$diff <- (body_size2$max) - (body_size2$min)
View(body_size2) # which have largest diff? Malurus_pulcherrimus, Charadrius_tricollaris, Thryothorus_ludovicianus

body_size_all <- gather(body_size, "dataset", "body_mass", c(2:12))
body_size_all <- body_size_all %>%
  dplyr::select(Species, dataset, body_mass)
# remove rows w NA brain volume
body_size_all2 <- body_size_all %>%
  filter(!is.na(body_mass))
body_size_all2 <- body_size_all2 %>%
  arrange(body_mass)

# species with biggest differences 

#Malurus_pulcherrimus 
which(body_size_all2$Species %in% "Malurus_pulcherrimus")
mp <- body_size_all2[c(2:307), ]
mp <- mp %>%
  dplyr::distinct(Species)
nrow(mp)-1  # minus itself, 140 species
((nrow(mp)-1)/nrow(body_size))*100 # 6.01%

#Charadrius_tricollaris
which(body_size_all2$Species %in% "Charadrius_tricollaris")
ac <- body_size_all2[c(1745:3393), ]
ac <- ac %>%
  dplyr::distinct(Species)
nrow(ac) # 759
((nrow(ac))/nrow(body_size))*100 # 32.56%

#Thryothorus_ludovicianus
which(body_size_all2$Species %in% "Thryothorus_ludovicianus")
tl <- body_size_all2[c(56:1170), ]
tl <- tl %>%
  dplyr::distinct(Species)
nrow(tl) # 429
((nrow(tl)-1)/nrow(body_size))*100 # 18.36
```

# figure 1
```{r echo = FALSE}
#tiff("fig_1.tiff", units="in", width=10, height=5, res=500)
panel <- grid.arrange(brain_size_plot,body_size_plot, nrow = 1) 
#dev.off()
```

# influence of variable inclusion, classification and source
```{r include = FALSE}
# subset to Sayol's data (most datapoints)
df <- corvides %>%
  filter(! is.na(LogBrainVol_SAYOL))
nrow(df) # 68 species

# remove 'EMAILED' (no reply when asking after variable info)
df <- df %>%
  filter(! SocialForaging_simp2 %in% "EMAILED")

# do results change if we treat partial migrants as resident?
df$Movement_partial_ed <- NA
for(i in 1:nrow(df)){
  if(df$Movement[i] == "Partial"){
    df$Movement_partial_ed[i] <- "Resident"} else {
      df$Movement_partial_ed[i] <- df$Movement[i]
    }
}

# do results change if we change facultative cooperative breeders or suspected but unknown cooperative breeders to independent breeders, rather than categorising them as cooperative breeders?

# Pyrrhocorax_graculus
# Corvus_brachyrhynchos
# Corvus_corone
# Artamus_cyanopterus
# Batis_capensis
# Bias_musicus
# Cinclosoma_cinnamomeum

df$Coop_ed <- NA
for(i in 1:nrow(df)){
  if(df$Species[i] %in% c("Pyrrhocorax_graculus","Corvus_brachyrhynchos","Corvus_corone","Artamus_cyanopterus","Batis_capensis","Bias_musicus","Cinclosoma_cinnamomeum")){
    df$Coop_ed[i] <- 0} else {
      df$Coop_ed[i] <- df$CooperativeBreeding_bin[i]
    }
}

## select variables to include
main <- df %>%
  dplyr::select(Species, LogBrainVol_SAYOL, LogBody_SAYOL, DietBreadth, TempVar, Movement, CooperativeBreeding_bin, SocialForaging_simp2, Coop_ed, Movement_partial_ed)
main$LogBrainVol_SAYOL <- as.numeric(main$LogBrainVol_SAYOL)
main$LogBody_SAYOL <- as.numeric(main$LogBody_SAYOL)
main$DietBreadth <- as.numeric(main$DietBreadth)
main$TempVar <- as.numeric(main$TempVar)
## create comparative object
main_comp<-comparative.data(tree2, main, names.col="Species", warn.dropped=T, force.root=T)
nrow(main_comp$data) # 59 species

# PPC dataframe
ppc <- df %>%
  dplyr::select(Species, LogBrainVol_SAYOL, LogBody_SAYOL, DietBreadth, PPC1, PPC2, Movement, CooperativeBreeding_bin, SocialForaging_simp2, Coop_ed, Movement_partial_ed)
## create comparative object
ppc$LogBrainVol_SAYOL <- as.numeric(ppc$LogBrainVol_SAYOL)
ppc$LogBody_SAYOL <- as.numeric(ppc$LogBody_SAYOL)
ppc$DietBreadth <- as.numeric(ppc$DietBreadth)
ppc$PPC1 <- as.numeric(ppc$PPC1)
ppc$PPC2 <- as.numeric(ppc$PPC2)
## create comparative object
ppc_comp<-comparative.data(tree2, ppc, names.col="Species", warn.dropped=T, force.root=T)
nrow(ppc_comp$data) # 46 species
```

## sih
```{r include = FALSE}
# main
sih <- pgls(LogBrainVol_SAYOL ~ LogBody_SAYOL + CooperativeBreeding_bin + SocialForaging_simp2, data=main_comp, lambda = "ML")

#model checks
lambda.profile<-pgls.profile(sih, "lambda")
plot(lambda.profile)
plot(sih)
sihvif <-  lm(LogBrainVol_SAYOL ~ LogBody_SAYOL + CooperativeBreeding_bin + SocialForaging_simp2, data=main_comp$data)
summary(sihvif)
vif(sihvif)

# summary
summary(sih)

# do results change if we treat partial migrants as resident?
sih_ed <- pgls(LogBrainVol_SAYOL ~ LogBody_SAYOL + Coop_ed + SocialForaging_simp2, data=main_comp, lambda = "ML")
summary(sih_ed)
summary(sih) # no change in results

# calculate CI
sih_ci <- sih_ed$sterr*1.96 
sih_upp_ci <- sih_ed$model$coef + sih_ci
sih_low_ci <- sih_ed$model$coef - sih_ci
```

## eih
```{r include = FALSE}
# main
eih <- pgls(LogBrainVol_SAYOL ~ LogBody_SAYOL + DietBreadth + TempVar + Movement, data=main_comp,lambda = "ML")

#model checks
lambda.profile<-pgls.profile(eih, "lambda")
plot(lambda.profile)
plot(eih)
eihvif <- lm(LogBrainVol_SAYOL ~ LogBody_SAYOL + DietBreadth + TempVar + Movement, data=main_comp$data)
summary(eihvif)
vif(eihvif)

# summary
summary(eih)

# do results change if we treat occasional cooperative breeders as non-cooperative breeders?
eih_ed <- pgls(LogBrainVol_SAYOL ~ LogBody_SAYOL + DietBreadth + TempVar + Movement_partial_ed, data=main_comp,lambda = "ML")
summary(eih_ed)
summary(eih)

# calculate CI
# eih
eih_ci <- eih$sterr*1.96 
eih_upp_ci <- eih$model$coef + eih_ci
eih_low_ci <- eih$model$coef - eih_ci
#eih_ed
eih_ed_ci <- eih_ed$sterr*1.96 
eih_ed_upp_ci <- eih_ed$model$coef + eih_ed_ci
eih_ed_low_ci <- eih_ed$model$coef - eih_ed_ci



# ppc
eih.ppc <- pgls(LogBrainVol_SAYOL ~ LogBody_SAYOL + DietBreadth + PPC1 + PPC2 + Movement, data=ppc_comp,lambda = "ML")

#model checks
lambda.profile<-pgls.profile(eih.ppc, "lambda")
plot(lambda.profile)
plot(eih.ppc)
eih.ppcvif <- lm(LogBrainVol_SAYOL ~ LogBody_SAYOL + DietBreadth + PPC1 + PPC2 + Movement, data=ppc_comp$data)
summary(eih.ppcvif)
vif(eih.ppcvif)

# summary
summary(eih.ppc)

# do results change when movement changed?
eih.ppc_ed <-  pgls(LogBrainVol_SAYOL ~ LogBody_SAYOL + DietBreadth + PPC1 + PPC2 + Movement_partial_ed, data=ppc_comp,lambda = "ML")
summary(eih.ppc_ed)

# calculate CIs
# eih
eih.ppc_ci <- eih.ppc$sterr*1.96 
eih.ppc_upp_ci <- eih.ppc$model$coef + eih.ppc_ci
eih.ppc_low_ci <- eih.ppc$model$coef - eih.ppc_ci
summary(eih.ppc)
#eih_ed
eih.ppc_ed_ci <- eih.ppc_ed$sterr*1.96 
eih.ppc_ed_upp_ci <- eih.ppc_ed$model$coef + eih.ppc_ed_ci
eih.ppc_ed_low_ci <- eih.ppc_ed$model$coef - eih.ppc_ed_ci
```

## combined
```{r include = FALSE}
# main
comb <- pgls(LogBrainVol_SAYOL ~ LogBody_SAYOL + DietBreadth + TempVar + Movement + CooperativeBreeding_bin + SocialForaging_simp2, data=main_comp,lambda = "ML")

#model checks
lambda.profile<-pgls.profile(comb, "lambda")
plot(lambda.profile)
plot(comb)
combvif <- lm(LogBrainVol_SAYOL ~ LogBody_SAYOL + DietBreadth + TempVar + Movement + CooperativeBreeding_bin + SocialForaging_simp2, data=main_comp$data)
summary(combvif)
vif(combvif)

#summary(comb)
summary(comb)

# do change in cooperative breeders change results?
# does change in movement change resultS?

comb_ed <- pgls(LogBrainVol_SAYOL ~ LogBody_SAYOL + DietBreadth + TempVar + Movement_partial_ed + Coop_ed + SocialForaging_simp2, data=main_comp,lambda = "ML")
summary(comb_ed)
summary(comb)

# calculate CIs
#comb
summary(comb)
comb_ci <- comb_ed$sterr*1.96 
comb_upp_ci <- comb$model$coef + comb_ci
comb_low_ci <- comb$model$coef - comb_ci

#comb_ed
summary(comb_ed)
comb_ed_ci <- comb_ed$sterr*1.96 
comb_ed_upp_ci <- comb_ed$model$coef + comb_ed_ci
comb_ed_low_ci <- comb_ed$model$coef - comb_ed_ci
```

# figure 2
```{r echo = FALSE}
main_comp$data$Movement <- as.character(main_comp$data$Movement)
main_comp$data$Movement[main_comp$data$Movement == "Migrant"] <- "zMigrant"
main_comp$data$Movement <- as.factor(main_comp$data$Movement)
main_comp$data$Movement <- as.character(main_comp$data$Movement)
main_comp$data$Movement[main_comp$data$Movement == "Partial"] <- "zPartial"
main_comp$data$Movement <- as.factor(main_comp$data$Movement)
ppc_comp$data$Movement <- as.character(ppc_comp$data$Movement)
ppc_comp$data$Movement[ppc_comp$data$Movement == "Migrant"] <- "zMigrant"
ppc_comp$data$Movement <- as.factor(ppc_comp$data$Movement)
ppc_comp$data$Movement <- as.character(ppc_comp$data$Movement)
ppc_comp$data$Movement[ppc_comp$data$Movement == "Partial"] <- "zPartial"
ppc_comp$data$Movement <- as.factor(ppc_comp$data$Movement)

#sih
est <- as.data.frame(sih$model$coef)
est$var <- rownames(est)
sih_plot <- dplyr::rename(est, sih = V1)

#eih
eih <- pgls(LogBrainVol_SAYOL ~ LogBody_SAYOL + DietBreadth + TempVar + Movement, data=main_comp,lambda = "ML")
summary(eih)
est <- as.data.frame(eih$model$coef)
est$var <- rownames(est)
eih_plot <- dplyr::rename(est, eih = V1)

#eih_ppc
eih.ppc <- pgls(LogBrainVol_SAYOL ~ LogBody_SAYOL + DietBreadth + PPC1 + PPC2 + Movement, data=ppc_comp,lambda = "ML")
summary(eih.ppc)
est <- as.data.frame(eih.ppc$model$coef)
est$var <- rownames(est)
eih.ppc_plot <- dplyr::rename(est, eih.ppc = V1)

#combined
comb <- pgls(LogBrainVol_SAYOL ~ LogBody_SAYOL + DietBreadth + TempVar + Movement + CooperativeBreeding_bin + SocialForaging_simp2, data=main_comp,lambda = "ML")
summary(comb)
est <- as.data.frame(comb$model$coef)
est$var <- rownames(est)
comb_plot <- dplyr::rename(est, comb = V1)

# combine all datasets
plot <- merge(sih_plot, eih_plot, by = "var", all.x = T, all.y = T)
plot <- merge(plot, eih.ppc_plot, by = "var", all.x = T, all.y = T)
plot <- merge(plot, comb_plot, by="var", all.x=T, all.y = T)

# change to long format
plot_long <- plot %>%
  pivot_longer(c("sih","eih","eih.ppc","comb") , names_to = "Model type", values_to = "estimate") # pivot_longer = new gather
plot_long
plot_long <- plot_long %>%
  dplyr::filter(! var %in% "(Intercept)") %>%
  dplyr::select(var, `Model type`, estimate)
plot_long$key <- paste0(plot_long$`Model type`, plot_long$var)

# add se

sihse <- as.data.frame(sih$sterr)
sihse$var <- rownames(sihse)
sihse <- dplyr::rename(sihse, sih = `sih$sterr`)

eihse <- as.data.frame(eih$sterr)
eihse$var <- rownames(eihse)
eihse <- dplyr::rename(eihse, eih = `eih$sterr`)

eih.ppcse <- as.data.frame(eih.ppc$sterr)
eih.ppcse$var <- rownames(eih.ppcse)
eih.ppcse <- dplyr::rename(eih.ppcse, eih.ppc = `eih.ppc$sterr`)

combse <- as.data.frame(comb$sterr)
combse$var <- rownames(combse)
combse <- dplyr::rename(combse, comb = `comb$sterr`)

plot <- merge(sihse, eihse, by = "var",all.x=T,all.y = T)
plot <- merge(plot, eih.ppcse, by = "var", all.x = T, all.y = T)
plot <- merge(plot, combse, by="var", all.x=T, all.y = T)

plot_longse <- plot %>%
  pivot_longer(c("sih","eih","eih.ppc","comb"), names_to = "Model type", values_to = "standard_error") # pivot_longer = new gather
plot_longse <- plot_longse %>%
  dplyr::filter(! var %in% "(Intercept)") %>%
  dplyr::select(var, `Model type`, standard_error)
plot_longse$key <- paste0(plot_longse$`Model type`, plot_longse$var)
plot_longse <- plot_longse %>%
  dplyr::select(standard_error,key)

plot_long3 <- merge(plot_long, plot_longse, by = "key", all.x=T, all.y = T)
names(plot_long3)

# dplyr::rename variables for plot
plot_long3[plot_long3 == "TempVar"] <- "Temperature Variation"
plot_long3[plot_long3 == "SocialForaging_simp2SG.nest"] <- "SF, nested small group"
plot_long3[plot_long3 == "SocialForaging_simp2SG"] <- "SF, non-nested small group"
plot_long3[plot_long3 == "SocialForaging_simp2S"] <- "SF, solitary"
plot_long3[plot_long3 == "SocialForaging_simp2P.nest"] <- "SF, nested pair"
plot_long3[plot_long3 == "MovementResident"] <- "Movement, resident"
plot_long3[plot_long3 == "MovementzPartial"] <- "Movement, partial"
plot_long3[plot_long3 == "MovementzMigrant"] <- "Movement, migrant"
plot_long3[plot_long3 == "MovementNomadic"] <- "Movement, nomadic"
plot_long3[plot_long3 == "LogBody_SAYOL"] <- "Body size (log)"
plot_long3[plot_long3 == "DietBreadth"] <- "Diet breadth"
plot_long3[plot_long3 == "CooperativeBreeding_bin"] <- "Cooperative breeding"

plot_long3[plot_long3 == "comb"] <- "combined"
plot_long3[plot_long3 == "eih.ppc"] <- "eih (ppc)"

# calculate CI
plot_long3$ci <- plot_long3$standard_error*1.96 
plot_long3$max_ci <- plot_long3$estimate + plot_long3$ci
plot_long3$min_ci <- plot_long3$estimate - plot_long3$ci

# colours
# The palette with grey:
palette <- c("#999999", "#E69F00", "#56B4E9", "#0072B2", "#D55E00", "#CC79A7")

#tiff("fig_2.tiff", units="in", width=6, height=5, res=500)
# plot
effect_plot <- ggplot(plot_long3, aes(estimate, var))+
  geom_point(aes(colour = `Model type`, shape = `Model type`), position=position_dodge(width=0.5)) +
  geom_errorbarh(mapping = aes(xmax = estimate + ci, xmin = estimate - ci, colour = `Model type`), position=position_dodge(width=0.5))+
  xlab("Estimate")+
  theme(legend.position = "bottom",legend.box="vertical",legend.margin=margin(),panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))+
  scale_colour_manual(values = palette)+
  scale_fill_manual(values = palette)+
  ylab("") +
  geom_vline(xintercept=0, linetype="dotted")+
  guides(color = guide_legend(nrow = 2, byrow = TRUE))
effect_plot
#dev.off()
```

# re-levelling models
## re-level: nested sociality (main)
```{r include = FALSE}
sih <- pgls(LogBrainVol_SAYOL ~ LogBody_SAYOL + CooperativeBreeding_bin + SocialForaging_simp2, data=main_comp, lambda = "ML")
summary(sih)

comb <- pgls(LogBrainVol_SAYOL ~ LogBody_SAYOL + DietBreadth + TempVar + Movement + CooperativeBreeding_bin + SocialForaging_simp2, data=main_comp,lambda = "ML")
summary(comb)

main_comp$data$SocialForaging_simp2 <- as.factor(main_comp$data$SocialForaging_simp2)
summary(main_comp$data$SocialForaging_simp2)
main_comp$data$SocialForaging_simp2 <- as.character(main_comp$data$SocialForaging_simp2)
main_comp$data$SocialForaging_simp2[main_comp$data$SocialForaging_simp2 == "P"] <- "zP"
main_comp$data$SocialForaging_simp2 <- as.factor(main_comp$data$SocialForaging_simp2)

sih.2 <- pgls(LogBrainVol_SAYOL ~ LogBody_SAYOL + CooperativeBreeding_bin + SocialForaging_simp2, data=main_comp, lambda = "ML")
summary(sih.2)

comb.2 <- pgls(LogBrainVol_SAYOL ~ LogBody_SAYOL + DietBreadth + TempVar + Movement + CooperativeBreeding_bin + SocialForaging_simp2, data=main_comp,lambda = "ML")
summary(comb.2)

summary(main_comp$SocialForaging_simp2)
main_comp$data$SocialForaging_simp2 <- as.character(main_comp$data$SocialForaging_simp2)
main_comp$data$SocialForaging_simp2[main_comp$data$SocialForaging_simp2 == "S"] <- "zS"
main_comp$data$SocialForaging_simp2 <- as.factor(main_comp$data$SocialForaging_simp2)

sih.3 <- pgls(LogBrainVol_SAYOL ~ LogBody_SAYOL + CooperativeBreeding_bin + SocialForaging_simp2, data=main_comp, lambda = "ML")
summary(sih.3)

comb.3 <- pgls(LogBrainVol_SAYOL ~ LogBody_SAYOL + DietBreadth + TempVar + Movement + CooperativeBreeding_bin + SocialForaging_simp2, data=main_comp,lambda = "ML")
summary(comb.3)

summary(main_comp$SocialForaging_simp2)
main_comp$data$SocialForaging_simp2 <- as.character(main_comp$data$SocialForaging_simp2)
main_comp$data$SocialForaging_simp2[main_comp$data$SocialForaging_simp2 == "SG"] <- "zSG"
main_comp$data$SocialForaging_simp2 <- as.factor(main_comp$data$SocialForaging_simp2)

sih.4 <- pgls(LogBrainVol_SAYOL ~ LogBody_SAYOL + CooperativeBreeding_bin + SocialForaging_simp2, data=main_comp, lambda = "ML")
summary(sih.4)

comb.4 <- pgls(LogBrainVol_SAYOL ~ LogBody_SAYOL + DietBreadth + TempVar + Movement + CooperativeBreeding_bin + SocialForaging_simp2, data=main_comp,lambda = "ML")
summary(comb.4)

```
```{r include = TRUE}
summary(sih)
summary(sih.2)
summary(sih.3)
summary(sih.4)

summary(comb)
summary(comb.2)
summary(comb.3)
summary(comb.4)

```

## re-level: movement (main)
```{r include = FALSE}
eih <- pgls(LogBrainVol_SAYOL ~ LogBody_SAYOL + DietBreadth + TempVar + Movement, data=main_comp,lambda = "ML")
summary(eih)
comb <- pgls(LogBrainVol_SAYOL ~ LogBody_SAYOL + DietBreadth + TempVar + Movement + CooperativeBreeding_bin + SocialForaging_simp2, data=main_comp,lambda = "ML")
summary(comb)

main_comp$data$Movement <- as.factor(main_comp$data$Movement) #4 levels
summary(main_comp$data$Movement)
main_comp$data$Movement <- as.character(main_comp$data$Movement)
main_comp$data$Movement[main_comp$data$Movement == "Migrant"] <- "zMigrant"
main_comp$data$Movement <- as.factor(main_comp$data$Movement)

eih.2 <- pgls(LogBrainVol_SAYOL ~ LogBody_SAYOL + DietBreadth + TempVar + Movement, data=main_comp,lambda = "ML")
summary(eih.2)

comb.2 <- pgls(LogBrainVol_SAYOL ~ LogBody_SAYOL + DietBreadth + TempVar + Movement + CooperativeBreeding_bin + SocialForaging_simp2, data=main_comp,lambda = "ML")
summary(comb.2)

summary(main_comp$Movement)
main_comp$data$Movement <- as.character(main_comp$data$Movement)
main_comp$data$Movement[main_comp$data$Movement == "Nomadic"] <- "zNomadic"
main_comp$data$Movement <- as.factor(main_comp$data$Movement)

eih.3 <- pgls(LogBrainVol_SAYOL ~ LogBody_SAYOL + DietBreadth + TempVar + Movement, data=main_comp,lambda = "ML")
summary(eih.3)

comb.3 <- pgls(LogBrainVol_SAYOL ~ LogBody_SAYOL + DietBreadth + TempVar + Movement + CooperativeBreeding_bin + SocialForaging_simp2, data=main_comp,lambda = "ML")
summary(comb.3)

summary(main_comp$Movement)
main_comp$data$Movement <- as.character(main_comp$data$Movement)
main_comp$data$Movement[main_comp$data$Movement == "Partial"] <- "zPartial"
main_comp$data$Movement <- as.factor(main_comp$data$Movement)

eih.4 <- pgls(LogBrainVol_SAYOL ~ LogBody_SAYOL + DietBreadth + TempVar + Movement, data=main_comp,lambda = "ML")
summary(eih.3)

comb.4 <- pgls(LogBrainVol_SAYOL ~ LogBody_SAYOL + DietBreadth + TempVar + Movement + CooperativeBreeding_bin + SocialForaging_simp2, data=main_comp,lambda = "ML")
summary(comb.3)

```
```{r include = TRUE}
summary(eih)
summary(eih.2)
summary(eih.3)
summary(eih.4)

summary(comb)
summary(comb.2)
summary(comb.3)
summary(comb.4)
```

## re-level: movement (ppc)
```{r include = FALSE}
eih <- pgls(LogBrainVol_SAYOL ~ LogBody_SAYOL + DietBreadth + PPC1 + PPC2 + Movement, data=ppc_comp,lambda = "ML")
summary(eih.2)

comb <- pgls(LogBrainVol_SAYOL ~ LogBody_SAYOL + DietBreadth + PPC1 + PPC2 + Movement + CooperativeBreeding_bin + SocialForaging_simp2, data=ppc_comp,lambda = "ML")
summary(comb.2)

ppc_comp$data$Movement <- as.factor(ppc_comp$data$Movement) #4 levels
summary(ppc_comp$data$Movement)
ppc_comp$data$Movement <- as.character(ppc_comp$data$Movement)
ppc_comp$data$Movement[ppc_comp$data$Movement == "Migrant"] <- "zMigrant"
ppc_comp$data$Movement <- as.factor(ppc_comp$data$Movement)

eih.2 <- pgls(LogBrainVol_SAYOL ~ LogBody_SAYOL + DietBreadth + PPC1 + PPC2 + Movement, data=ppc_comp,lambda = "ML")
summary(eih.2)

comb.2 <- pgls(LogBrainVol_SAYOL ~ LogBody_SAYOL + DietBreadth + PPC1 + PPC2 + Movement + CooperativeBreeding_bin + SocialForaging_simp2, data=ppc_comp,lambda = "ML")
summary(comb.2)

summary(ppc_comp$Movement)
ppc_comp$data$Movement <- as.character(ppc_comp$data$Movement)
ppc_comp$data$Movement[ppc_comp$data$Movement == "Nomadic"] <- "zNomadic"
ppc_comp$data$Movement <- as.factor(ppc_comp$data$Movement)

eih.3 <- pgls(LogBrainVol_SAYOL ~ LogBody_SAYOL + DietBreadth + PPC1 + PPC2 + Movement, data=ppc_comp,lambda = "ML")
summary(eih.3)

comb.3 <- pgls(LogBrainVol_SAYOL ~ LogBody_SAYOL + DietBreadth + PPC1 + PPC2 + Movement + CooperativeBreeding_bin + SocialForaging_simp2, data=ppc_comp,lambda = "ML")
summary(comb.3)

summary(ppc_comp$Movement)
ppc_comp$data$Movement <- as.character(ppc_comp$data$Movement)
ppc_comp$data$Movement[ppc_comp$data$Movement == "Partial"] <- "zPartial"
ppc_comp$data$Movement <- as.factor(ppc_comp$data$Movement)

eih.4 <- pgls(LogBrainVol_SAYOL ~ LogBody_SAYOL + DietBreadth + PPC1 + PPC2 + Movement, data=ppc_comp,lambda = "ML")
summary(eih.4)

comb.4 <- pgls(LogBrainVol_SAYOL ~ LogBody_SAYOL + DietBreadth + PPC1 + PPC2 + Movement + CooperativeBreeding_bin + SocialForaging_simp2, data=ppc_comp,lambda = "ML")
summary(comb.4)


```
```{r include = TRUE}
summary(eih)
summary(eih.2)
summary(eih.3)
summary(eih.4)

summary(comb)
summary(comb.2)
summary(comb.3)
summary(comb.4)
```