---
title: "Analysis code - Smith & Garcia 2021"
author: "Morphological responses to food limitation in a habitat forming sea urchin, Strongylocentrotus purpuratus"
date: "3/16/2021"
output: html_document:
  toc:yes
  pdf_document:
    highlight: tango
    toc: yes
    toc_depth: 6
---

**Required data**
#This R markdown file requires the data file 'Smith_data_file.csv' available on Dryad at XXXX

**Required Packages**

```{r load_libraries, message=FALSE}
library(ggplot2)
library(ggsignif)
library(dplyr)
library(Hmisc)
library(ggpubr)
library(gridExtra)
library(grid)
library(ggpmisc)
library(betareg)
```

**Load Data for ANCOVAs and shape plots**
```{r load_data}
#set appropriate file directory
urchin.data<-data.frame(read.csv(file.choose()))
```

**Data Legend for urchin.data**
Year = survey year
Site.Number = survey site
Density = mean urchin density (per m^2) across all 16 1x1m quadrats within a given site
Brown = mean % cover of brown algae from 16 1x1m quadrats
Encrusting = mean % cover of encrusting red and coralline algae from 16 1x1m quadrats
Red = mean % cover of foliose red algae from 16 1x1m quadrats
TH = mean urchin test height (mm)
TD = mean urchin test diameter (mm)
LL = mean urchin lantern length (mm)
LW = mean urchin lantern width (mm)
LI = mean urchin lantern index
GI = mean urchin gonad index
Zone = Barren (BAR), or Forest (FOR)
GWM = mean gonad wet mass (g)
Log_GWM = log transformed GWM
Log_TD = log transformed TD

**Data Legend for zone.data**
see meta data above. Field listed in zone.data are summarized as mean values for each test diameter and zone. 


**Part 1 - plot shape data**
```{r}
ll.lw<-ggplot(urchin.data, aes(x = LW, y = LL, shape=ZONE, colour = ZONE, fill=ZONE)) +  
  geom_point() +
  geom_smooth(method='lm') +
  scale_fill_manual(values=c("#7C5295", "#228B22"))+
  scale_colour_manual(values=c("#7C5295", "#228B22"))+
  theme(# Hide panel borders and remove grid lines
    panel.background = element_blank(),
    panel.border = element_blank(),
    panel.grid.major = element_line(colour = c("white"), size = c(0.33, 0.2)),
    panel.grid.minor = element_blank(),
    # Change axis line
    axis.line = element_line(colour = "black"),
    axis.text.x=element_text(size=16, color=c("black", NA)),
    plot.title = element_text(hjust=0.5, face="bold"),
    axis.text.y=element_text(size=16, color = c("black",NA)),
    text=element_text(size=16),
    legend.position="none"
  ) +
  scale_y_continuous(breaks = seq(8, 14, by=0.5),lim=c(8, 14))+
  scale_x_continuous(breaks = seq(7, 13, by=0.5),lim=c(7, 13))+
  xlab("lantern width (mm)")+
  ylab("lantern length (mm)")+
  stat_regline_equation(
    aes(label =  paste(..eq.label.., ..adj.rr.label.., sep = "~~~~")),
    formula = formula,
    size=6
  )

ll.lw



td.th<-ggplot(urchin.data, aes(x = TD, y = TH, shape=ZONE, colour = ZONE, fill=ZONE)) +  
  geom_point() +
  geom_smooth(method='lm') +
  scale_fill_manual(values=c("#7C5295", "#228B22"))+
  scale_colour_manual(values=c("#7C5295", "#228B22"))+
  theme(# Hide panel borders and remove grid lines
    panel.background = element_blank(),
    panel.border = element_blank(),
    panel.grid.major = element_line(colour = c("white"), size = c(0.33, 0.2)),
    panel.grid.minor = element_blank(),
    # Change axis line
    axis.line = element_line(colour = "black"),
    axis.text.x=element_text(size=16, color=c("black", NA)),
    plot.title = element_text(hjust=0.5, face="bold"),
    axis.text.y=element_text(size=16, color = c("black",NA)),
    text=element_text(size=16),
    legend.position="none"
  ) +
  scale_y_continuous(breaks = seq(13, 21, by=0.5),lim=c(13, 21))+
  scale_x_continuous(breaks = seq(30, 45, by=1),lim=c(30, 45))+
  xlab("test diameter (mm)")+
  ylab("test height (mm)")+
  stat_regline_equation(
    aes(label =  paste(..eq.label.., ..adj.rr.label.., sep = "~~~~")),
    formula = formula,
    size=6
  )

td.th

figure.2 <- ggarrange(
  ll.lw, td.th,
  labels = "AUTO",
  ncol = 2, nrow = 1)

figure.2
```

**Part 2 - variation in sea urchin morphometrics**
#The first test here is an ANCOVA used for evaluating variation in lantern length across habitat types (barren, forest) with test diamter as a covariate

```{r}
#Run ANCOVA for lantern length
ll.ANCOVA <- lm(urchin.data$LL~urchin.data$TD+relevel(urchin.data$ZONE, ref="FOR"), data=urchin.data)
summary(ll.ANCOVA)
```

#The second test is an ANCOVA used for evaluating variation in gonad mass across habitat types (barren, forest) with test diamter as a covariate

```{r}
#Gonad index ANCOVA 
GWM.ANCOVA <- lm(urchin.data$GWM~urchin.data$TD+relevel(urchin.data$ZONE, ref="FOR"), data=urchin.data)
summary(GWM.ANCOVA)
```


#plot relationship between gonad mass, lantern length, and test diameter
**Load zone summarized data for plotting**
```{r load_data}
#set appropriate file directory
zone.data<-data.frame(read.csv(file.choose()))
```

```{r}
formula = y ~ x


ll.plot<-ggplot(zone.data, aes(x = Log_TD, y = Log_LL, shape=ZONE, colour = ZONE, fill=ZONE)) +  
  geom_point() +
  geom_smooth(method = "lm") +
  scale_fill_manual(values=c("#7C5295", "#228B22"))+
  scale_colour_manual(values=c("#7C5295", "#228B22"))+
  theme(# Hide panel borders and remove grid lines
    panel.background = element_blank(),
    panel.border = element_blank(),
    panel.grid.major = element_line(colour = c("white"), size = c(0.33, 0.2)),
    panel.grid.minor = element_blank(),
    # Change axis line
    axis.line = element_line(colour = "black"),
    axis.text.x=element_text(size=14, color=c("black", NA)),
    plot.title = element_text(hjust=0.5, face="bold"),
    axis.text.y=element_text(size=14, color = c(NA, "black")),
    text=element_text(size=14),
    legend.position="none"
  ) +
  scale_x_continuous()+
  scale_y_continuous(breaks = seq(2.2,3, by=0.2), lim=c(2.2,3))+
  xlab("log(test diameter [mm])")+
  ylab("log(lantern length [mm])")+
  stat_regline_equation(
    aes(label =  paste(..eq.label.., ..adj.rr.label.., sep = "~~~~")),
    formula = formula
  )
ll.plot

```

#plot relationship between gonad mass and test diameter
```{r}
GWM.plot<-ggplot(zone.data, aes(x = Log_TD, y = Log_GWM, shape = ZONE, colour = ZONE, fill=ZONE)) +  
  geom_smooth(method = "lm") +
  geom_point()+
  #stat_smooth(method='lm', aes(ymin = ifelse(..ymin.. < 0, 0, ..ymin..)), 
              #alpha = .3) +
  
  scale_fill_manual(values=c("#7C5295", "#228B22"),labels=c("barren","forest"))+
  scale_colour_manual(values=c("#7C5295", "#228B22"),labels=c("barren","forest"))+
  theme(# Hide panel borders and remove grid lines
    panel.background = element_blank(),
    panel.border = element_blank(),
    panel.grid.major = element_line(colour = c("white"), size = c(0.33, 0.2)),
    panel.grid.minor = element_blank(),
    # Change axis line
    axis.line = element_line(colour = "black"),
    axis.text.x=element_text(size=14, color=c("black", NA)),
    plot.title = element_text(hjust=0.5, face="bold"),
    axis.text.y=element_text(size=14, color = c(NA, "black")),
    text=element_text(size=14),
    legend.position= "none", #edit full plot legend here and extract using function below
    legend.title=element_blank(),
    legend.text = element_text(size=18)
  ) +
  scale_x_continuous()+
  scale_y_continuous(breaks = seq(-2,3, by=1), limits =c(-2,3))+
  xlab("log(test diameter [mm])")+
  ylab("log(gonad wet weight [g])")+
  stat_regline_equation(
    aes(label =  paste(..eq.label.., ..adj.rr.label.., sep = "~~~~")),
    formula = formula
  ) 
GWM.plot



figure.3 <- ggarrange(
                    ll.plot, 
                    GWM.plot,
                    labels = "AUTO",
                    ncol = 1, nrow = 2,
                    heights=c(4,4),
                    align = "v")
          
figure.3

```

**Part 3 - continuous models using lantern index and gonad index as a function of algae % cover**

#Construct linear models and plots
```{r}
LI.brown<-ggplot(urchin.data, aes(x = Brown, y = LI)) +  
  geom_point() +
  geom_smooth(method = "lm", color="black", fill = "black") +
  theme(# Hide panel borders and remove grid lines
    panel.background = element_blank(),
    panel.border = element_blank(),
    panel.grid.major = element_line(colour = c("white"), size = c(0.33, 0.2)),
    panel.grid.minor = element_blank(),
    # Change axis line
    axis.line = element_line(colour = "black"),
    axis.text.x=element_text(size=14, color=c("black", NA)),
    plot.title = element_text(hjust=0.5, face="bold"),
    axis.text.y=element_text(size=14, color = c(NA,"black")),
    text=element_text(size=14),
    legend.position="none"
  ) +
  scale_y_continuous(breaks = seq(0.26, 0.37, by=0.01),lim=c(0.26, 0.37))+
  scale_x_continuous(breaks = seq(0, 40, by=5),lim=c(0, 40))+
  xlab("")+
  ylab("")+
  stat_poly_eq(formula = y ~ x,
               aes(label = paste(..eq.label.., ..rr.label.., ..p.value.label.., sep = "*','~")),
               parse=TRUE,
               label.x.npc = "left",
               vstep = 0.05,
               size=4) 
LI.brown

LI.encrusting<-ggplot(urchin.data, aes(x = Encrusting, y = LI)) +  
  geom_point() +
  geom_smooth(method = "lm", color="black", fill = "black") +
  theme(# Hide panel borders and remove grid lines
    panel.background = element_blank(),
    panel.border = element_blank(),
    panel.grid.major = element_line(colour = c("white"), size = c(0.33, 0.2)),
    panel.grid.minor = element_blank(),
    # Change axis line
    axis.line = element_line(colour = "black"),
    axis.text.x=element_text(size=14, color=c("black", NA)),
    plot.title = element_text(hjust=0.5, face="bold"),
    axis.text.y=element_text(size=14, color = c(NA,"black")),
    text=element_text(size=14),
    legend.position="none"
  ) +
  scale_y_continuous(breaks = seq(0.26, 0.37, by=0.01),lim=c(0.26, 0.37))+
  scale_x_continuous(breaks = seq(0, 100, by=5),lim=c(0, 100))+
  xlab("")+
  ylab("")+
  stat_poly_eq(formula = y ~ x,
               aes(label = paste(..eq.label.., ..rr.label.., ..p.value.label.., sep = "*','~")),
               parse=TRUE,
               label.x.npc = "left",
               vstep = 0.05,
               size=4) 

LI.encrusting


LI.Red<-ggplot(urchin.data, aes(x = Red, y = LI)) +  
  geom_point() +
  geom_smooth(method = "lm", color="black", fill = "black") +
  theme(# Hide panel borders and remove grid lines
    panel.background = element_blank(),
    panel.border = element_blank(),
    panel.grid.major = element_line(colour = c("white"), size = c(0.33, 0.2)),
    panel.grid.minor = element_blank(),
    # Change axis line
    axis.line = element_line(colour = "black"),
    axis.text.x=element_text(size=14, color=c("black", NA)),
    plot.title = element_text(hjust=0.5, face="bold"),
    axis.text.y=element_text(size=14, color = c(NA,"black")),
    text=element_text(size=14),
    legend.position="none"
  ) +
  scale_y_continuous(breaks = seq(0.26, 0.37, by=0.01),lim=c(0.26, 0.37))+
  scale_x_continuous(breaks = seq(0, 100, by=5),lim=c(0, 100))+
  xlab("")+
  ylab("")+
  stat_poly_eq(formula = y ~ x,
               aes(label = paste(..eq.label.., ..rr.label.., ..p.value.label.., sep = "*','~")),
               parse=TRUE,
               label.x.npc = "left",
               vstep = 0.05,
               size=4) 

LI.Red



GI.brown<-ggplot(urchin.data, aes(x = Brown, y = GI)) +  
  geom_point() +
  geom_smooth(method = "lm", color="black", fill = "black") +
  theme(# Hide panel borders and remove grid lines
    panel.background = element_blank(),
    panel.border = element_blank(),
    panel.grid.major = element_line(colour = c("white"), size = c(0.33, 0.2)),
    panel.grid.minor = element_blank(),
    # Change axis line
    axis.line = element_line(colour = "black"),
    axis.text.x=element_text(size=14, color=c("black", NA)),
    plot.title = element_text(hjust=0.5, face="bold"),
    axis.text.y=element_text(size=14, color = c(NA,"black")),
    text=element_text(size=14),
    legend.position="none"
  ) +
  scale_y_continuous(breaks = seq(0, 16, by=1),lim=c(0, 16))+
  scale_x_continuous(breaks = seq(0, 40, by=5),lim=c(0, 40))+
  xlab("")+
  ylab("")+
  stat_poly_eq(formula = y ~ x,
               aes(label = paste(..eq.label.., ..rr.label.., ..p.value.label.., sep = "*','~")),
               parse=TRUE,
               label.x.npc = "left",
               vstep = 0.05,
               size=4) 


GI.brown

GI.encrusting<-ggplot(urchin.data, aes(x = Encrusting, y = GI)) +  
  geom_point() +
  stat_smooth(method='lm', aes(ymin = ifelse(..ymin.. < 0, 0, ..ymin..)), 
              alpha = .3, color='black', fill='black') +
  theme(# Hide panel borders and remove grid lines
    panel.background = element_blank(),
    panel.border = element_blank(),
    panel.grid.major = element_line(colour = c("white"), size = c(0.33, 0.2)),
    panel.grid.minor = element_blank(),
    # Change axis line
    axis.line = element_line(colour = "black"),
    axis.text.x=element_text(size=14, color=c("black",NA)),
    plot.title = element_text(hjust=0.5, face="bold"),
    axis.text.y=element_text(size=14, color = c("black",NA)),
    text=element_text(size=14),
    legend.position="none"
  ) +
  scale_y_continuous(breaks = seq(0, 16, by=1),lim=c(0, 16))+
  scale_x_continuous(breaks = seq(0, 100, by=5),lim=c(0, 100))+
  xlab("")+
  ylab("")+
  stat_poly_eq(formula = y ~ x,
               aes(label = paste(..eq.label.., ..rr.label.., ..p.value.label.., sep = "*','~")),
               parse=TRUE,
               label.x.npc = "left",
               vstep = 0.05,
               size=4) 

GI.encrusting


GI.Red<-ggplot(urchin.data, aes(x = Red, y = GI)) +  
  geom_point() +
  geom_smooth(method = "lm", color="black", fill = "black") +
  theme(# Hide panel borders and remove grid lines
    panel.background = element_blank(),
    panel.border = element_blank(),
    panel.grid.major = element_line(colour = c("white"), size = c(0.33, 0.2)),
    panel.grid.minor = element_blank(),
    # Change axis line
    axis.line = element_line(colour = "black"),
    axis.text.x=element_text(size=14, color=c("black", NA)),
    plot.title = element_text(hjust=0.5, face="bold"),
    axis.text.y=element_text(size=14, color = c(NA,"black")),
    text=element_text(size=14),
    legend.position="none"
  ) +
  scale_y_continuous(breaks = seq(0, 16, by=1),lim=c(0, 16))+
  scale_x_continuous(breaks = seq(0, 100, by=5),lim=c(0, 100))+
  xlab("")+
  ylab("")+
  stat_poly_eq(formula = y ~ x,
                aes(label = paste(..eq.label.., ..rr.label.., ..p.value.label.., sep = "*','~")),
                parse=TRUE,
                label.x.npc = "left",
                vstep = 0.05,
                size=4) 


GI.Red






figure.4 <- ggarrange(
                    LI.encrusting, GI.encrusting,
                    LI.brown,GI.brown,
                    LI.Red, GI.Red,
                    labels = "AUTO",
                    ncol = 2, nrow = 3)

figure.4

annotate_figure(figure,
                bottom = text_grob("cover of algae from photoquadrats (%)", size=18, vjust=-0.9),
                left = text_grob("lantern index", size = 18, rot = 90, hjust=0, vjust=1),
                right = text_grob("gonad index (%)", size = 18, rot=-90, hjust=0.8, vjust=1))
```

**Part 4 - beta regression**
#evaluate algae cover as a function of urchin density
```{r}
y.transform.betareg <- function(y){
    n.obs <- sum(!is.na(y))
    (y*(n.obs-1) + 0.5) / n.obs
}

#scale percent cover to proportation (0:1 response required for betareg)
Brown = (urchin.data$Brown)/100
Red = (urchin.data$Red)/100
Encrusting = (urchin.data$Encrusting)/100

gy_logit_brown <- betareg(y.transform.betareg(Brown)~urchin.data$Density)
gy_logit_red <- betareg(y.transform.betareg(Red)~urchin.data$Density)
gy_logit_enc <- betareg(y.transform.betareg(Encrusting)~urchin.data$Density)
gy_logit_LI <- betareg(y.transform.betareg(urchin.data$LI)~urchin.data$Density)



#summary reports
summary(gy_logit_brown)
summary(gy_logit_red)
summary(gy_logit_enc)
```

#plot results from betareg
```{r}
algae.plot<-ggplot(urchin.data, aes(x = urchin.data$Density, y = urchin.data$Brown)) +
         #geom_point(size=4, shape=21) +
         scale_fill_grey() +
         geom_line(aes(y = predict(gy_logit_brown, urchin.data), 
                       colour = "Brown"), size=1.2) +
         geom_line(aes(y = predict(gy_logit_red, urchin.data), 
                       colour = "Red"), size=1.2) +
         geom_line(aes(y = predict(gy_logit_enc, urchin.data), 
                       colour = "Encrusting"), size=1.2) +
         scale_colour_manual(values=c(Brown = "brown", Red = "red", Encrusting = "purple"))+
         scale_linetype_manual("", values = c("solid"))+
         ylab(expression(Algae~percent~cover))+
         xlab(expression(Mean~purple~urchin~density~(per~m^2)))
         theme_bw()

print(algae.plot + theme_bw(base_size = 20) + theme(panel.grid.minor = element_blank(), legend.title=element_blank()) + ylab(expression(Proportion~of~cover~algae)) + xlab(expression(Mean~purple~urchin~density~(per~m^2))))
```


**End Report**