---
title: "Field Study"
author: "Amanda Franklin, Alberto Rivera, Justin Robbins, Jan Pechenik"
date: "12/12/2020"
output: pdf_document
---

## Does BMI vary through winter or with environmental conditions?

This script runs the field study analysis for the article "Body Mass Index does 
not decline during winter for the sedentary marine gastropod Crepidula 
fornicata" by Amanda Franklin, Alberto Rivera, Justin Robbins and Jan Pechenik.

## Set up

Load required packages

```{r setup}
# Packages

library(tidyverse)
library(lubridate)
library(car)
library(lme4)
library(plyr)
library(lsmeans)
library(lmtest)

```

Load data and process data

```{r data}
# read in data

dat <- read.csv("../dat/data_BMI.csv")
sstchl <- read.csv("../dat/data_envvar.csv")

# Data clean/add variables
dat2 <- dat %>%
  dplyr::rename(weight = Dry.tissue.weight.mg, length = Shell.Length.mm) %>%
  mutate(date = dmy(Date), 
         year = year(date),
         month = month(date),
         lengthcen = length-min(length, na.rm=T),
         lengthcm = length * 0.1,
         dayyear = yday(date),
         totalmonth = ((year - 2009)*12) + month,
         fyear = as.factor(year),
         BMI = weight/(lengthcm^3),
         sex = case_when(
           length <= 25 ~ "Male",
           length >= 30 ~ "Female",
           length = TRUE ~ "Ambiguous"
         )) %>%
  filter(!is.na(length)) 

datall <- merge(dat2, sstchl, by = "Date")
```

Look at data and create plot S1.

```{r data exploration}
# Sample size:
length(datall$BMI)
tapply(datall$BMI, datall$sex, FUN = length)

## Have a look at variance:
boxplot(datall$BMI) 

## Check normality/spread
ggplot(datall,aes(x=BMI)) + 
  geom_histogram() +
  theme_bw()   # seems good.

## Plot all individuals

plot1 <- ggplot(dat = datall, aes(length, weight, colour = sex)) + 
  geom_point(size = 0.8) +
  theme_classic() +
  scale_colour_manual(values = c("grey50", "gold1", "tomato1")) +
  xlab("Shell length (mm)") + ylab("Dry tissue weight (mg)") +
  theme(axis.title.x = element_text(face="bold", 
                                    size=13),  #change axes size/bold
        axis.text.x  = element_text(size=11, colour="black"),
        axis.title.y = element_text(face="bold", size = 13, vjust=1),
        axis.text.y = element_text (size=11, colour="black"),
        axis.line.x = element_line(color="black", size = 0.5),
        axis.line.y = element_line(color="black", size = 0.5),
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        legend.title=element_blank(), 
        legend.text = element_text(size = 11),
        legend.justification=c(0,0), legend.position=c(0.01,0.65)) 
plot1


ggsave(plot1, file = "../fig/lengthweight.jpg", width = 100, height = 80,
       units = "mm", dpi=300)

```

## PCA with environmental variables
Look at the environmental variables and run PCA due to high correlation between
variables

```{r analysis PCA}

## Look at data
# Relationship with month:
plot(datall$SST_Monthly ~ datall$month)
plot(datall$Chlorophyll_monthly ~ datall$month)
plot(datall$AirT_MonthlyAvg ~ datall$month)

# Correlation between variables
plot(datall$SST_Monthly ~ datall$Chlorophyll_monthly) # moderate to high correlation
plot(datall$SST_Monthly ~ datall$AirT_MonthlyAvg) # high correlation
plot(datall$AirT_MonthlyAvg ~ datall$Chlorophyll_monthly) # moderate to high correlation

## PCA
# Due to high correlation, run PCA and use PC as predictor
# Scale variables:
datall$SST_scale <- scale(datall$SST_Monthly)[,1]
datall$Chl_scale <- scale(datall$Chlorophyll_monthly)[,1]
datall$AirT_scale <- scale(datall$AirT_MonthlyAvg)[,1]

# PCA:
pc <- prcomp(~ SST_scale + Chl_scale + AirT_scale, datall)
# look at loadings
pc
# look at variance
summary(pc)
# Look at eigenvalues
screeplot(pc, type="lines")
# add to dataset
datall$pc1 <- pc$x[,1]

## Plot PCA by month
plot_pcsmth <- ggplot(datall, aes(x=month, y = pc1)) +
  geom_point(size=1) + 
  scale_y_continuous(name = "PC1") + #,
  theme_classic() + xlab("Month") +
  theme(
    axis.title.y = element_text(size=13)) +
  theme(axis.title.x = element_text(face="bold", 
                                    size=13),  #change axes size/bold
        axis.text.x  = element_text(size=11, colour="black", angle = 45,
                                    vjust = 0.5),
        axis.title.y = element_text(face="bold", size = 13, vjust=1),
        axis.text.y = element_text (size=11, colour="black"),
        axis.line.x = element_line(color="black", size = 0.5),
        axis.line.y = element_line(color="black", size = 0.5)) +
  scale_x_continuous(breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9 ,10, 11, 12), 
                     label = c("Jan", "Feb", "Mar", "Apr", "May", "Jun", 
                     "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"))
plot_pcsmth

ggsave(plot_pcsmth, file = "../fig/PCsMonth.jpg", width = 100, height = 80,
       units = "mm", dpi=1000)


```
## Analysis: BMI and Environmental Variables

Does BMI correlate with PC1 (temperature and chlorophyll a)? Run a linear model
and use likelihood ratio tests to determine significance of each variable.

```{r analysis PCA BMI}

m0 <- lm(BMI ~ pc1*sex + totalmonth, datall)
m1 <- lm(BMI ~ pc1 + sex + totalmonth, datall)
m2 <- lm(BMI ~ pc1 + sex, datall)
m3 <- lm(BMI ~ pc1 + totalmonth, datall)
m4 <- lm(BMI ~ sex + totalmonth, datall)

plot(m0)
summary(m0)

lrtest(m0, m1)
lrtest(m1, m2)
lrtest(m1, m3)
lrtest(m1, m4)



plot_pc1 <- ggplot(datall, aes(x = pc1, y = BMI, colour = sex)) +
  geom_jitter(width = 0.02, alpha = 0.3) +
  scale_colour_manual(values = c("grey50", "gold1", "tomato1")) +
  stat_function(fun = function(x)(-0.558187 *(x) + 9.319730), 
                geom = "line", col = "gold3",
                xlim = c(-1.940793, 3.671914), size = 1) +
  stat_function(fun = function(x)(-0.558187 *(x) + 9.319730  + 0.796702), 
                geom = "line", col = "tomato3",
                xlim = c(-1.940793, 3.671914), size = 1) +
  theme_classic() + xlab("PC1") +
  theme(axis.title.x = element_text(face="bold", 
                                    size=13),  #change axes size/bold
        axis.text.x  = element_text(size=11, colour="black"),
        axis.title.y = element_text(face="bold", size = 13, vjust=1),
        axis.text.y = element_text (size=11, colour="black"),
        axis.line.x = element_line(color="black", size = 0.5),
        axis.line.y = element_line(color="black", size = 0.5),
        legend.title=element_blank(), 
        legend.text = element_text(size = 11),
        legend.justification=c(0,0), legend.position=c(0.6,0.7))
plot_pc1

ggsave(plot_pc1, file = "../fig/BMIPC1.jpg", width = 100, height = 80,
       units = "mm", dpi=1000)


```


```{r winter data}

# Filter out the data we want and add a variable indicating which winter
# data is from

dat_winter <- datall %>%
  filter((month <= 4 & year %in% c(2010, 2017, 2018)) | 
        (month >= 10 & year %in% c(2009, 2016, 2017)),
        sex != "N/A") %>%
  mutate(winteryear = case_when(
    year %in% c(2009, 2010) ~ "2009/10",
    year == 2016 ~ "2016/17",
    year == 2018 ~ "2017/18",
    year == 2017 & month <= 4 ~ "2016/17",
    year == 2017 & month >= 10 ~ "2017/18"
  ),
        month2 = case_when(
          month <= 4 ~ 12 + month, 
          month >= 10 ~ month
        ))

```

```{r change through winter data exploration}

## Sample size
tapply(dat_winter$BMI, dat_winter$month2, FUN = length)

## Have a look at variance:
boxplot(dat_winter$BMI ~ dat_winter$month2) # looks fine

# Check:
var1 <- dat_winter %>%
  filter(month2 == 15)
var2 <- dat_winter %>%
  filter(month2 == 16)

var((var1$BMI))/var((var2$BMI)) #not bad (much less than 4)


## Check normality/spread

ggplot(dat_winter, aes(x=BMI)) + 
  geom_histogram() +
  facet_grid(~month2)+theme_bw()   # seems good.

```

```{r winter weight loss analysis}
m2_all <- lm(BMI ~ month2 * sex + winteryear, dat = dat_winter)

# Check model fit
plot(m2_all)  

# Results
summary(m2_all)

m2_all1 <- lm(BMI ~ sex + winteryear + month2, dat = dat_winter)
m2_all2 <- lm(BMI ~ winteryear + month2, dat = dat_winter)
m2_all3 <- lm(BMI ~ sex + winteryear, dat = dat_winter)
m2_all4 <- lm(BMI ~ sex + month2, dat = dat_winter)

lrtest(m2_all, m2_all1)
lrtest(m2_all1, m2_all2)
lrtest(m2_all1, m2_all3)
lrtest(m2_all1, m2_all4)

summary(m2_all)



```

```{r plot change through winter}
summarydat <- dat_winter %>%
  group_by(winteryear, month2) %>%
  dplyr::summarise(mean = mean(BMI),
            sd = sd(BMI),
            n = length(BMI), 
            se = sd/sqrt(n),
            BMI = mean) 


plot_winter <- ggplot(dat_winter, aes(x=month2, y = BMI, 
                                group = winteryear, colour = winteryear)) + 
  geom_point(position = position_jitterdodge(), alpha = 0.6, size = 0.9) +
  geom_errorbar(dat = summarydat, 
                aes(x = month2, 
                    ymin = BMI - se, 
                    ymax = BMI + se, 
                    group = winteryear),
                colour = "black",  position = position_dodge(width = 0.9), width = 0.2) +
  geom_point(dat = summarydat, aes(x=month2, y = BMI, 
                                   group = winteryear, fill = winteryear), 
             size = 2, colour = "black", shape = 21, 
             position = position_dodge(width = 0.9)) +
  scale_color_manual(values = c("deepskyblue4", "springgreen4", "goldenrod3"), 
                     labels = c("2009-2010", "2016-2017", "2017-2018")) +
  scale_fill_manual(values = c("deepskyblue4", "springgreen4", "goldenrod3"),
                    labels = c("2009-2010", "2016-2017", "2017-2018")) +
  scale_x_continuous(breaks = c(10, 11, 12, 13, 14, 15, 16), 
                     label = c("Oct", "Nov", "Dec", "Jan", "Feb", "Mar", "Apr")) +
  theme_classic() +
  theme(legend.title = element_blank()) +
  xlab("Month") +
  theme(axis.title.x = element_text(face="bold", 
                                    size=13),  #change axes size/bold
        axis.text.x  = element_text(size=11, colour="black"),
        axis.title.y = element_text(face="bold", size = 13, vjust=1),
        axis.text.y = element_text (size=11, colour="black"),
        axis.line.x = element_line(color="black", size = 0.5),
        axis.line.y = element_line(color="black", size = 0.5),
        legend.title=element_blank(), 
        legend.text = element_text(size = 11),
        legend.justification=c(0,0), legend.position=c(0.4,0.7),
        strip.text.x = element_text(size=13, hjust = 0),
        legend.background = element_blank()) +
    guides(color = guide_legend(override.aes = list(size = 2) ) )
plot_winter


ggsave(plot_winter, file = "../fig/WinterWeight_2022.jpg", width = 100, height = 80,
       units = "mm", dpi=1000)

```
