---
title: "Code for modelling hierarchal variation in metabolic scaling exponents"
author: "Fonti Kar"
date: "10/12/2018"
output: html_document
---

```{r setup, include=FALSE}
rm(list=ls())
sessionInfo()
#R version 3.4.2 (2017-09-28)
#Platform: x86_64-apple-darwin15.6.0 (64-bit)
#Running under: macOS High Sierra 10.13.6

packages.need <- c("ggplot2", "tidyverse", "dplyr", "magrittr", "MCMCglmm", "brms", "RColorBrewer")
lapply(packages.need, require, character.only = TRUE)
```

## Read in long format data first

```{r load data}
data <- read.csv("R/RMDs/Long_data_for_analysis.csv")
str(data)
dim(data)
```

#Calculating among and within-individual mass deviations
For more details: van de Pol, M., & Wright, J. (2009). A simple method for distinguishing within- versus between-subject effects using mixed models. Animal Behaviour, 77(3), 753???758. http://doi.org/10.1016/j.anbehav.2008.11.006

```{r}
ms_data <- data %>% 
  select(id, samp_period, incb_temp, log.co2pmin, log.mass) %>%
  group_by(id) %>% 
  mutate(AMIlogmass = mean(log.mass, na.rm = T)) %>%  #Among ID effect
  mutate(WIlogmass = round(log.mass - AMIlogmass, 2)) %>%  #Within ID effect (Deviation of an individual's mass relative to its own mean)
  arrange(id, samp_period, incb_temp) %>% 
  as.data.frame() 
str(ms_data)
```

#Code to modelling the scaling relationship between metabolic rate and among- and within- individual mass effects and to test temeprature dependence
```{r}
#With interaction with temperature
mod.ms.1a <- brm(log.co2pmin ~ as.factor(incb_temp)*AMIlogmass + as.factor(incb_temp)*WIlogmass + (1 + WIlogmass|id), 
            data = ms_data,
            family = gaussian(),
            cores = 3)

mod.ms.1a <- add_ic(mod.ms.1a, ic = c("loo","waic"))
#saveRDS(mod.ms.1a, "R/RMDs/output/mod.ms.1a") 

#No interaction with temperature
mod.ms.1b <- brm(log.co2pmin ~ as.factor(incb_temp) + AMIlogmass + WIlogmass + (1 + WIlogmass|id), 
            data = ms_data,
            family = gaussian(),
            cores = 3)

mod.ms.1b <- add_ic(mod.ms.1b, ic = c("loo","waic"))
#saveRDS(mod.ms.1b, "R/RMDs/output/mod.ms.1b") 

#Read in model output
mod.ms.1a <- readRDS("R/RMDs/output/mod.ms.1a")
mod.ms.1b <- readRDS("R/RMDs/output/mod.ms.1b")

#Diagnostics
plot(mod.ms.1a)
summary(mod.ms.1a)

plot(mod.ms.1b)
summary(mod.ms.1b)

#Compare wAIC values
waic(mod.ms.1a, mod.ms.1b)   

#Compare wloo values
loo(mod.ms.1a, mod.ms.1b)  
```

### Code to generate Figure 5A and Figure S4
```{r}
#Trying with entire posterior
mod.ms.1.fixed <- as.matrix(posterior_samples(mod.ms.1a, "^b_")) #Selecting the fixed effects

#Creating table of mass-scaling exponents for all temperatures at the AmongID and WithinID level
hier.ms.dat <- as.data.frame(matrix(ncol = 5, nrow = 6*2))
colnames(hier.ms.dat) <- c("Temp", "Type", "Est", "Lower", "Upper")
hier.ms.dat$Temp <- rep(seq(22,32, by =2), 2)
hier.ms.dat$Type <- rep(c("Among ID", "Within ID"), each = 6)

#Among ID effects
#22C
hier.ms.dat[1,3:5] <- posterior_summary(mod.ms.1.fixed)[7,c(1,3:4)]

#24C
hier.ms.dat[2,3] <- mean(mod.ms.1.fixed[,7] + mod.ms.1.fixed[,9])
hier.ms.dat[2,4:5] <-HPDinterval(as.mcmc(mod.ms.1.fixed[,7] + mod.ms.1.fixed[,9]))

#26C
hier.ms.dat[3,3] <- mean(mod.ms.1.fixed[,7] + mod.ms.1.fixed[,10])
hier.ms.dat[3,4:5] <- HPDinterval(as.mcmc(mod.ms.1.fixed[,7] + mod.ms.1.fixed[,10]))

#28
hier.ms.dat[4,3] <- mean(mod.ms.1.fixed[,7] + mod.ms.1.fixed[,11])
hier.ms.dat[4,4:5] <-  HPDinterval(as.mcmc(mod.ms.1.fixed[,7] + mod.ms.1.fixed[,11]))

#30
hier.ms.dat[5,3] <- mean(mod.ms.1.fixed[,7] + mod.ms.1.fixed[,12])
hier.ms.dat[5,4:5] <- HPDinterval(as.mcmc(mod.ms.1.fixed[,7] + mod.ms.1.fixed[,12]))

#32
hier.ms.dat[6,3] <- mean(mod.ms.1.fixed[,7] + mod.ms.1.fixed[,13])
hier.ms.dat[6,4:5] <- HPDinterval(as.mcmc(mod.ms.1.fixed[,7] + mod.ms.1.fixed[,13]))

#Within ID effects
#22C
#22C
hier.ms.dat[7,3] <- mean(mod.ms.1.fixed[,8])
hier.ms.dat[7,4:5] <- HPDinterval(as.mcmc(mod.ms.1.fixed[,8]))

#24C
hier.ms.dat[8,3] <- mean(mod.ms.1.fixed[,8] + mod.ms.1.fixed[,14])
hier.ms.dat[8,4:5] <- HPDinterval(as.mcmc(mod.ms.1.fixed[,8] + mod.ms.1.fixed[,14]))

#26C
hier.ms.dat[9,3] <- mean(mod.ms.1.fixed[,8] + mod.ms.1.fixed[,15])
hier.ms.dat[9,4:5] <- HPDinterval(as.mcmc(mod.ms.1.fixed[,8] + mod.ms.1.fixed[,15]))

#28
hier.ms.dat[10,3] <- mean(mod.ms.1.fixed[,8] + mod.ms.1.fixed[,16])
hier.ms.dat[10,4:5] <- HPDinterval(as.mcmc(mod.ms.1.fixed[,8] + mod.ms.1.fixed[,16]))

#30
hier.ms.dat[11,3] <- mean(mod.ms.1.fixed[,8] + mod.ms.1.fixed[,17])
hier.ms.dat[11,4:5] <- HPDinterval(as.mcmc(mod.ms.1.fixed[,8] + mod.ms.1.fixed[,17]))

#32
hier.ms.dat[12,3] <- mean(mod.ms.1.fixed[,8] + mod.ms.1.fixed[,18])
hier.ms.dat[12,4:5] <- HPDinterval(as.mcmc(mod.ms.1.fixed[,8] + mod.ms.1.fixed[,18]))

#Theme for all ggplots
my_theme <- theme_bw() + 
  theme(legend.position = "none", 
        legend.text = element_text(size = 10),
        legend.title = element_blank(),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.text=element_text(size=14),
        axis.title=element_text(size=18),
        strip.text =element_text(size=18)) 

filter(hier.ms.dat, Type == "Among ID") %>%  ggplot(aes(y = Est, x = Temp)) +
  geom_errorbar(aes(colour = Type, ymin = Lower, ymax = Upper), width = 0, position = position_dodge(1.5)) +
  geom_point(aes(fill = Type, shape = Type), size = 4, position = position_dodge(1.5)) + 
  geom_hline(aes(yintercept = 0.83), linetype = 2) +
  scale_color_manual(values = c(rep("black",2))) + 
  scale_shape_manual(values= c(24, 23)) + 
  scale_fill_manual(values= c("black", "grey")) +
  labs(x = "Temperature", y = "Mass-scaling exponent estimate") + 
  theme_bw() +
  scale_x_continuous(breaks = seq(22,32, by = 2)) +
  scale_y_continuous(breaks = seq(-3.5,3.5, by = 0.5)) +
  my_theme +
  theme(legend.position = c(0.90, 0.10),
        legend.box.background = element_rect(colour = "black")) 
 
#Figure S4
 ggplot(data = hier.ms.dat, aes(y = Est, x = Temp)) +
  geom_errorbar(aes(colour = Type, ymin = Lower, ymax = Upper), width = 0, position = position_dodge(1.5)) +
  geom_point(aes(fill = Type, shape = Type), size = 4, position = position_dodge(1.5)) + 
  geom_hline(aes(yintercept = 0.83), linetype = 2) +
  scale_color_manual(values = c(rep("black",2))) + 
  scale_shape_manual(values= c(24, 23)) + 
  scale_fill_manual(values= c("black", "grey")) +
  labs(x = "Temperature", y = "Mass-scaling exponent estimate") + 
  theme_bw() +
  scale_x_continuous(breaks = seq(22,32, by = 2)) +
  scale_y_continuous(breaks = seq(-3.5,3.5, by = 0.5)) +
  my_theme +
  theme(legend.position = c(0.90, 0.10),
        legend.box.background = element_rect(colour = "black")) 
 
#Fig. 5B 
#Choosing unique colours
n <- 42
qual_col_pals = brewer.pal.info[brewer.pal.info$category == 'qual',]
col_vector = unlist(mapply(brewer.pal, qual_col_pals$maxcolors, rownames(qual_col_pals)))

set.seed(10)
my_cols <- sample(col_vector, n, replace = F) ; pie(rep(1,n), col=my_cols)
 
set.seed(100)
mysubset <- sample(data$id, 20, replace = F) 

data %>% filter(id %in% mysubset) %>%
  ggplot(aes(x = log.mass  , y = log.co2pmin, group = "id")) + 
  geom_point(shape = 21, aes(fill = id), size = 1.5, color = "black") + 
  geom_line(aes(group = id), stat="smooth", method = "lm", alpha = 0.3, lwd = 0.5) +
  geom_line(stat = "smooth", method = "lm", alpha = 1, lwd = 1) +
  facet_grid(~ incb_temp) + 
  scale_fill_manual(values= my_cols) +
  scale_x_continuous(breaks = c(-0.20, 0, 0.20, 0.40)) +
  labs(x = expression(Log~Mass~(g)), y = expression(Log~Metabolic~rate~(VCO[2]~min^{-1}))) + 
  my_theme
```

###Code to generate Table 4
```{r}
Table4 <- data.frame(matrix(nrow = 6, ncol = 10))
colnames(Table4) <- c("Temperature", rep(c("Estimate", "Lower", "Upper"),3))
Table4$Temperature <- seq(22,32, by = 2)

#AmongID
Table4[,2:4] <- hier.ms.dat[1:6,3:5]

#WithinID
Table4[,5:7] <- hier.ms.dat[7:12,3:5]

#Difference
Table4[,8:10] <- hier.ms.dat[13:18,3:5]

Table4[,c(2:10)] <- round(Table4[,c(2:10)],2)

write.csv(Table4, row.names = F, "R/RMDs/output/tabs/Table4.csv")

```
