---
title: "Figure 1 — Hardison 2022"
author: "Emily Hardison"
date: "10/5/2022"
output: html_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
```{r}
############
# reading in packages
############
library(tidyverse)
library(lemon)
library(here)
#Linear Models
library(lme4) # fitting non-linear model
#library(lmerTest) # get P-values from non-linear Model
############
# adding function for calculating standard error
############
se<- function(x, na.rm=FALSE) {
if (na.rm) x <- na.omit(x)
sqrt(var(x)/length(x))
}
############
# reading in csv
############
acc_cap <- read_csv(here::here("hardison2022_abt_TPC_data.csv")) %>%
mutate(diet_treatment = as.factor(diet_treatment))%>%
mutate(fish_id = as.factor(fish_id))%>%
mutate(timepoint = as.factor(timepoint)) %>%
mutate(diet_fullname = case_when(diet_treatment == "Herb"~"Herbivorous",
diet_treatment == "Carn"~"Carnivorous",
diet_treatment == "Omni"~"Omnivorous"))
```
Acclimation capacity model fits
```{r}
############
# Fitting mixed effect models using lmer package for BIC comparison
############
#Linear:
m.linear.1 = lmer(as.numeric(fh_bpm)~acute_temperature + (1|fish_id)+ (1|replicate_number),
data=acc_cap, REML=FALSE)
#Linear:
m.linear.2 = lmer(as.numeric(fh_bpm)~acute_temperature + timepoint + (1|fish_id)+ (1|replicate_number),
data=acc_cap, REML=FALSE)
#Linear:
m.linear.3 = lmer(as.numeric(fh_bpm)~acute_temperature + as.factor(diet_treatment) + (1|fish_id)+ (1|replicate_number),
data=acc_cap, REML=FALSE)
#Linear:
m.linear.4 = lmer(as.numeric(fh_bpm)~acute_temperature + timepoint + as.factor(diet_treatment) + (1|fish_id)+ (1|replicate_number),
data=acc_cap, REML=FALSE) # REML means is restricting parameter space of data for the maximum likelihood estimate of your model. Essentially it limits model, and we opted to exclude it.
#Linear w/ interaction:
model.linear.int1 = lmer(as.numeric(fh_bpm)~acute_temperature + timepoint*as.factor(diet_treatment) + (1|fish_id)+ (1|replicate_number),
data=acc_cap, REML=FALSE)
#Linear w/ interaction:
model.linear.int2 = lmer(as.numeric(fh_bpm)~acute_temperature*timepoint*as.factor(diet_treatment) + (1|fish_id)+ (1|replicate_number),
data=acc_cap, REML=FALSE)
#non-linear, only acute acute_temperature:
m0 = lmer(as.numeric(fh_bpm)~poly(acute_temperature, 3) + (1|fish_id)+ (1|replicate_number),
data=acc_cap, REML=FALSE)
#3rd order poly:
m1 = lmer(as.numeric(fh_bpm)~poly(acute_temperature,3) + timepoint + as.factor(diet_treatment) + (1|fish_id)+ (1|replicate_number),
data=acc_cap, REML=FALSE)
#3rd order poly:
m2 = lmer(as.numeric(fh_bpm)~poly(acute_temperature,3) * timepoint * as.factor(diet_treatment) + (1|fish_id)+ (1|replicate_number),
data=acc_cap, REML=FALSE)
#3rd order poly:
m3 = lmer(as.numeric(fh_bpm)~poly(acute_temperature,3) * timepoint + (1|fish_id)+ (1|replicate_number),
data=acc_cap, REML=FALSE)
#3rd order poly:
m4 = lmer(as.numeric(fh_bpm)~poly(acute_temperature,3) * as.factor(diet_treatment) + (1|fish_id)+ (1|replicate_number),
data=acc_cap, REML=FALSE)
#3rd order poly:
m5 = lmer(as.numeric(fh_bpm)~poly(acute_temperature,3) * timepoint + as.factor(diet_treatment) + (1|fish_id)+ (1|replicate_number),
data=acc_cap, REML=FALSE)
#4th order poly:
m6 = lmer(as.numeric(fh_bpm)~poly(acute_temperature,4) * timepoint + as.factor(diet_treatment) + (1|fish_id)+ (1|replicate_number),
data=acc_cap, REML=FALSE)
#4th order poly:
m7 = lmer(as.numeric(fh_bpm)~poly(acute_temperature,4) * timepoint * as.factor(diet_treatment) + (1|fish_id)+ (1|replicate_number),
data=acc_cap, REML=FALSE)
#4th order poly:
m8 = lmer(as.numeric(fh_bpm)~poly(acute_temperature,4) + timepoint + as.factor(diet_treatment) + (1|fish_id)+ (1|replicate_number),
data=acc_cap, REML=FALSE)
#4th order poly:
m9 = lmer(as.numeric(fh_bpm)~poly(acute_temperature,4) * timepoint + (1|fish_id)+ (1|replicate_number),
data=acc_cap, REML=FALSE)
#4th order poly:
m10 = lmer(as.numeric(fh_bpm)~poly(acute_temperature,4) + timepoint * as.factor(diet_treatment) + (1|fish_id)+ (1|replicate_number),
data=acc_cap, REML=FALSE)
############
# Generating table of BIC comparisons
############
bic_table_all <- BIC(m0, m1, m2, m3, m4, m5,
m6, m7, m8, m9, m10,
m.linear.1,
m.linear.2,
m.linear.3,
m.linear.4,
model.linear.int1,
model.linear.int2) %>%
rownames_to_column(var = "model.name")
bic_table_all <- bic_table_all[order(bic_table_all$BIC), ]
bic_table_all$delta <- abs(bic_table_all$BIC[1] - bic_table_all$BIC)
bic_table_all <- bic_table_all%>%
mutate(delta = round(delta, 3))
##############
# The best fit model was m2
##############
```
Generating curves for plot of acclimation capacity
```{r}
############
# Looking at residual plot for abt data from model-- comparing normality of residuals between linear and model including diet, temp effects
############
# Linear model
plot(m.linear.1)
hist(resid(m.linear.1))
# Best fit Model
plot(m2)
hist(resid(m2))
############
# Fitting curves to data to generate Figure 1
############
# Fitting thermal performance curves (TPC) to the data for fish held two weeks at 12°C. Only fitting data to the average T[ARR] + 1xSD to not fit the curve past where the data exists
# Make dummy dataframe and fill with model info
predict.data_t12<- as.data.frame(expand.grid(timepoint = "0",
diet_treatment = levels(acc_cap$diet_treatment),
acute_temperature = seq(11.5, 30.7, by = 0.75))) # did curve predict by 0.75 intervals so that I could make linetype = acclimation temperature in plot.
predict.data_t12 <- predict.data_t12%>%
mutate(
predict_bpm = predict(m2, newdata = predict.data_t12, re.form = NA),# Adding in predicted model values to dummy df
treatment = paste(diet_treatment, timepoint, sep = "_"))
# Fitting Thermal performance curve (TPC) to the data for fish held two weeks at 20°C. Only fitting data to the average T[ARR] + 1xSD to not fit the curve past where the data exists
# Make dummy dataframe and fill with model info
predict.data_t20<- as.data.frame(expand.grid(timepoint = c("14"), diet_treatment = levels(acc_cap$diet_treatment), acute_temperature = seq(17.3, 34, by = 0.75))) # did curve predict by 0.75 intervals so that I could make linetype = acclimation temperature. when the interval <0.5 then the line always looks solid because of how close the values are to one another
predict.data_t20<-predict.data_t20 %>%
mutate(
predict_bpm = predict(m2, newdata = predict.data_t20, re.form = NA),# Adding in predicted model values
treatment = paste(diet_treatment, timepoint, sep = "_"))
# Bind data from different temperatures to make sure all data is there and the different temps are accounted for
predict.data <- rbind(predict.data_t12, predict.data_t20)
############
# Generating Figure 1
############
figure1 <- ggplot(data = acc_cap, aes(x = acute_temperature, y = fh_bpm, color = diet_treatment, fill = diet_treatment), na.rm = T)+
geom_line(aes(group = fish_id), alpha = 0.25)+
geom_line(data = predict.data, aes(x = acute_temperature, y = predict_bpm,
group = treatment,
color = diet_treatment,
linetype = timepoint), na.rm =T, size = 2)+ # Adding in Model -- full interaction
geom_point(size = 1.35, aes(shape = timepoint, alpha = timepoint))+
scale_y_continuous(limits = c(0, 260), expand = c(0,0), breaks = c(50, 100, 150, 200, 250))+
scale_x_continuous(limits = c(10, 38), expand = c(0,0), breaks = c(14, 18, 22, 26, 30, 34))+
labs(y = expression(italic(f)[Hmax]~(bpm)),
x = expression(paste("Temperature (", degree, "C)")),
color = "Treatment") +
scale_fill_manual(values = c("#BA4E24", "#368B73", "#D9AB39"))+
scale_color_manual(values = c("#BA4E24", "#368B73", "#D9AB39"))+
theme_classic()+
theme(axis.text = element_text(size=10, color = "black"),
axis.title.y= element_text(size=12, color = "black"),
axis.title.x= element_text(size=12),
legend.position = "none",
axis.line = element_line(colour = 'black', size = 0.5),
strip.background = element_rect(
color="white", fill="white"),
strip.text = element_blank())+
scale_shape_manual(values = c(19, 1))+
scale_alpha_manual(values = c(0.35, 0.55))+
facet_rep_grid(~diet_treatment)+
geom_text(mapping = aes(x = 15,
y = 245,
label = diet_fullname,
size = 1))
############
# Saving Figure 1
############
ggsave(figure1, filename = here::here( "figure1_Hardison2022.png"), width = 10, height = 4)
#ggsave(figure1, filename = here::here( "figure1_Hardison2022.pdf"), width = 10, height = 4)
```