R Session Info

sessionInfo()
## R version 3.4.3 (2017-11-30)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.3
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## loaded via a namespace (and not attached):
##  [1] compiler_3.4.3  backports_1.1.2 magrittr_1.5    rprojroot_1.3-2
##  [5] tools_3.4.3     htmltools_0.3.6 yaml_2.1.16     Rcpp_0.12.15   
##  [9] stringi_1.1.6   rmarkdown_1.8   knitr_1.19      stringr_1.2.0  
## [13] digest_0.6.15   evaluate_0.10.1

EPSCoR SCTC Aquatic Ecology

Juvenile salmon diet data - 2015 - Analysis in progress

Draft in progress 2/7/18. Quick graphic exploration of 2015 diet data.

#require packages
library(googlesheets)
library(tidyverse)
library(forcats)
suppressMessages(library(dplyr))
# detach("package:MASS", unload=TRUE) #pkg masks some dplyr functions

# Clear environment
rm(list=ls())
###################################################################
# Create Stacked Bar Charts of Diet Proportions by Sampling Event #
#
# a.) Finer-level catgories
# b.) coarser-level categories
###################################################################


# identify gs sheet object
SCTC2015 <- gs_title("2015 EPSCoR SCTC.xlsx")

#read in worksheet containing diet data
diet <- SCTC2015 %>% 
  gs_read(ws = "2015 Diet Contents Data",col_names=TRUE) %>%
  
#rename needed columns  
rename(prcsdt=Process_Date,
      sample.id=Sample.ID,
      spp=Fish_Species,
      river=River,
      reach=Reach,
      site=Site,
      date=Date,
      sampledate=Sample_Date,
      taxkey=Taxon_Keyed,
      prey=Prey_Type_Used,
      preycat=PreyCategory,
      dm_tot=Total_Prey_Dry_Mass_mg,
      sample.event=Sample.Event,
      event.num=sample.event.num,
      TA=TA_Final) %>%

#select desired columns
select(sample.id,sample.event,river,
       reach,site,sampledate,preycat,spp,dm_tot,TA,event.num) 

#filter out non-diet samples and blank data
diet <- diet %>% filter(site!="DR1",
       !is.na(spp),
       !is.na(sample.event),
       !is.na(preycat),
       !is.na(dm_tot),
       !is.na(event.num)
       )
    
#transform columns to desired classes
diet <- transform(diet, 
                  dm_tot = as.numeric(dm_tot),
                  sample.event = as.factor(sample.event),
                  event.num = as.factor(event.num),
                  spp = as.factor(spp))  %>%
  mutate(season = fct_recode(event.num,         #code sampling events as seasons
                          "Early Summer"="1", 
                          "Mid-Summer" = "2",
                          "Late Summer" ="3"))

#create sum of pooled prey DM for each prey category by fish species and sample event 
dm_sum_by_event <- diet %>%
  group_by(river,reach,sample.event,spp,preycat,TA,season) %>%
  summarise(sum(dm_tot)) %>%
  rename(dm_sum=`sum(dm_tot)`) 
  
#create table of summed dry mass for each fish species by sampling event
dm_sum_spp <- dm_sum_by_event %>%
  group_by(sample.event,spp) %>%
  summarise(pool_sum = sum(dm_sum)) 

#join pooled prey DM values to table grouped by fish species and prey category
dm_sum_final <- dm_sum_spp %>%
  left_join(dm_sum_by_event,by=c("sample.event","spp")) %>%
  mutate(pct = (dm_sum/pool_sum)*100) %>%
  filter(sample.event!="#N/A") 

#create table of quantity of Chinook and coho non-empty diet samples from each SITE VISIT
diet_sample_count <- SCTC2015 %>% 
  gs_read(ws = "C 2015 Diet & LW Data",col_names=TRUE) %>%
  filter(DS=="Y") %>%
  mutate(season = fct_recode(as.factor(sample.event.num), 
                             "Early Summer"="1", 
                             "Mid-Summer" = "2",
                             "Late Summer" ="3")) %>%
  group_by(sample.event,season,spp) %>%
  summarise(diet_count=n())

#create table of quantity of Chinook and coho non-empty diet samples from each SEASON
diet_sample_season <- SCTC2015 %>% 
  gs_read(ws = "C 2015 Diet & LW Data",col_names=TRUE) %>%
  filter(DS=="Y") %>%
  mutate(season = fct_recode(as.factor(sample.event.num), 
                             "Early Summer"="1", 
                             "Mid-Summer" = "2",
                             "Late Summer" ="3")) %>%
  group_by(river,season,spp) %>%
  summarise(diet_count=n()) %>%
  mutate(event_label = paste(as.character(spp), "( n=", as.character(diet_count),")"))


#join diet sample dry mass calculations by sampling event with fish count table by event
#e.g., how many fish of each species were diet-sampled at each sampling event?

#data is lumped by sampling event
dm_sum_final_by_event <- dm_sum_final %>%
  left_join(diet_sample_count, 
            by=c("season","spp", "sample.event")) %>% 
  mutate(event_label = paste("n =", as.character(diet_count)))

#join diet sample dry mass calculations by sampling event with fish count table by season
#e.g., how many fish of each species were diet-sampled in each season?

#data is lumped for each drainage by `season`
dm_sum_final_by_season <- dm_sum_final_by_event %>%
  group_by(season, river, spp, preycat, pct) %>%
  summarise(dry_mass_pooled = sum(dm_sum)) 



# create objects that contain only one species per drainage, summarizing 2015
# diet data

### Coho ###

# results for Beaver Creek 2015 only, Coho only
dm_sum_final_beaver_coho_2015 <- dm_sum_final_by_event %>%
  filter(river=="Beaver Creek", spp=="Coho") 
  
# results for Russian River 2016 only, Coho only
dm_sum_final_russian_coho_2015 <- dm_sum_final_by_event %>%
  filter(river=="Russian River", spp=="Coho") 

# results for Ptarmigan Creek 2015 only, Coho only
dm_sum_final_ptarmigan_coho_2015 <- dm_sum_final_by_event %>%
  filter(river=="Ptarmigan Creek", spp=="Coho")

# results for Kenai River 2015 only, Coho only
dm_sum_final_kenai_coho_2015 <- dm_sum_final_by_event %>%
  filter(river=="Kenai River", spp=="Coho")

### Chinook ###

# results for Beaver Creek 2015 only, Chinook only
dm_sum_final_beaver_chinook_2015 <- dm_sum_final_by_event %>%
  filter(river=="Beaver Creek", spp=="Chinook") 

# results for Russian River 2016 only, Chinook only
dm_sum_final_russian_chinook_2015 <- dm_sum_final_by_event %>%
  filter(river=="Russian River", spp=="Chinook") 

# results for Ptarmigan Creek 2015 only, Chinook only
dm_sum_final_ptarmigan_chinook_2015 <- dm_sum_final_by_event %>%
  filter(river=="Ptarmigan Creek", spp=="Chinook")

# results for Kenai River 2015 only, Chinook only
dm_sum_final_kenai_chinook_2015 <- dm_sum_final_by_event %>%
  filter(river=="Kenai River", spp=="Chinook")


##################
#                #
# Visualizations #
#                #
##################





#######################################
# Chinook/Coho diet proportions 2015 
#
# a.) FINE-level diet data (~5 prey categories)
#
# pooled by SITE and SEASON
#
########################################

# Create format for displaying fine-level diet data by site/season 
finelevel15 <- ggplot() +
  labs(x="",y="Percent Dry Mass (%)") +
  scale_fill_discrete(name="Prey Category",
                      breaks=c("Fish",
                               "InvertAquatic",
                               "InvertTerrestrial_AqOrigin",
                               "InvertTerrestrial_TerOrigin",
                               "InvertTerrestrial_UnkOrigin",
                               "InvertUnknown",
                               "FishEggs"),
                      labels=c("Fish",
                               "Aquatic Invertebrate",
                               "Terrestrial Invertebrate,\nAquatic Origin",
                               "Terrestrial Invertebrate,\nTerrestrial Origin       ",
                               "Terrestrial Invertebrate,\nUnknown Origin",
                               "Unknown Invertebrate",
                               "Fish Eggs")) +
  theme(legend.key.size = unit(1.7, 'lines')) +
  theme(plot.title = element_text(size = 24, face = "bold"),
        axis.text.x  = element_text(size=14),
        axis.title.y = element_text(size = 14, face = "bold")) +
  facet_grid(season~reach) +
  theme(strip.text.x = element_text(size=10, face="bold"), 
        strip.text.y = element_text(size=10, face="bold"))


################
# Beaver Creek #
################

#### Coho ####

#Create faceted bar charts of dry mass proportions by species and sampling event
finelevel15 + geom_bar(aes(y=pct, x=spp, fill=preycat), 
                       data = dm_sum_final_beaver_coho_2015,
                       stat="identity") +
  geom_text(data=dm_sum_final_beaver_coho_2015, 
            aes(x=1, y=4, label=event_label),
            vjust = 0.5) +
  ggtitle("Beaver Creek 2015\n Juvenile Coho Diet (% Dry Mass)")

#### Chinook ####

#Create faceted bar charts of dry mass proportions by species and sampling event
finelevel15 + geom_bar(aes(y=pct, x=spp, fill=preycat), 
                       data = dm_sum_final_beaver_chinook_2015,
                       stat="identity") +
  geom_text(data=dm_sum_final_beaver_chinook_2015, 
            aes(x=1, y=4, label=event_label),
            vjust = 0.5) +
  ggtitle("Beaver Creek 2015\n Juvenile Chinook Diet (% Dry Mass)")

#################
# Russian River #
#################

#### Coho ####

#Create faceted bar charts of dry mass proportions by species and sampling event
finelevel15 + geom_bar(aes(y=pct, x=spp, fill=preycat), 
                       data = dm_sum_final_russian_coho_2015,
                       stat="identity") +
  geom_text(data=dm_sum_final_russian_coho_2015, 
            aes(x=1, y=4, label=event_label),
            vjust = 0.5) +
  ggtitle("Russian River 2015\n Juvenile Coho Diet (% Dry Mass)")

#### Chinook ####

#Create faceted bar charts of dry mass proportions by species and sampling event
finelevel15 + geom_bar(aes(y=pct, x=spp, fill=preycat), 
                       data = dm_sum_final_russian_chinook_2015,
                       stat="identity") +
  geom_text(data=dm_sum_final_russian_chinook_2015, 
            aes(x=1, y=4, label=event_label),
            vjust = 0.5) +
  ggtitle("Russian River 2015\n Juvenile Chinook Diet (% Dry Mass)")

###################
# Ptarmigan Creek #
###################

#### Coho ####

#Create faceted bar charts of dry mass proportions by species and sampling event
finelevel15 + geom_bar(aes(y=pct, x=spp, fill=preycat), 
                       data = dm_sum_final_ptarmigan_coho_2015,
                       stat="identity") +
  geom_text(data=dm_sum_final_ptarmigan_coho_2015, 
            aes(x=1, y=4, label=event_label),
            vjust = 0.5) +
  ggtitle("Ptarmigan Creek 2015\n Juvenile Coho Diet (% Dry Mass)")

#### Chinook ####

#Create faceted bar charts of dry mass proportions by species and sampling event
finelevel15 + geom_bar(aes(y=pct, x=spp, fill=preycat), 
                       data = dm_sum_final_ptarmigan_chinook_2015,
                       stat="identity") +
  geom_text(data=dm_sum_final_ptarmigan_chinook_2015, 
            aes(x=1, y=4, label=event_label),
            vjust = 0.5) +
  ggtitle("Ptarmigan Creek 2015\n Juvenile Chinook Diet (% Dry Mass)")

###############
# Kenai River #
###############

#### Coho ####

#Create faceted bar charts of dry mass proportions by species and sampling event
finelevel15 + geom_bar(aes(y=pct, x=spp, fill=preycat), 
                       data = dm_sum_final_kenai_coho_2015,
                       stat="identity") +
  geom_text(data=dm_sum_final_kenai_coho_2015, 
            aes(x=1, y=4, label=event_label),
            vjust = 0.5) +
  ggtitle("Kenai River 2015\n Juvenile Coho Diet (% Dry Mass)")

#### Chinook ####

#Create faceted bar charts of dry mass proportions by species and sampling event
finelevel15 + geom_bar(aes(y=pct, x=spp, fill=preycat), 
                       data = dm_sum_final_kenai_chinook_2015,
                       stat="identity") +
  geom_text(data=dm_sum_final_kenai_chinook_2015, 
            aes(x=1, y=4, label=event_label),
            vjust = 0.5) +
  ggtitle("Kenai River 2015\n Juvenile Chinook Diet (% Dry Mass)")

 # ------------------------------------------------------------------------ #


#######################################
#
# b.) COARSE-level prey category data
# 
# e.g., Terrestrial vs. Aquatic source
#
#######################################

# Create format for displaying diet data by site/season 
coarselevel15 <- ggplot() +
  labs(x="",y="Percent Dry Mass (%)") +
  scale_fill_discrete(name="Prey Category",
                      breaks=c("T","A","M","UNK"),
                      labels=c("Terrestrial",
                               "Aquatic",
                               "Marine",
                               "Unknown")) +
  theme(legend.key.size = unit(1.7, 'lines')) +
  theme(plot.title = element_text(size = 24, face = "bold"),
        axis.text.x  = element_text(size=14),
        axis.title.y = element_text(size = 14, face = "bold")) +
  facet_grid(season~reach) +
  theme(strip.text.x = element_text(size=10, face="bold"), 
        strip.text.y = element_text(size=10, face="bold"))


################
# Beaver Creek #
################

#### Coho ####

#Create faceted bar charts of dry mass proportions by species and sampling event
coarselevel15 + geom_bar(aes(y=pct, x=spp, fill=TA, order=TA), 
                       data = dm_sum_final_beaver_coho_2015,
                       stat="identity") +
  geom_text(data=dm_sum_final_beaver_coho_2015, 
            aes(x=1, y=4, label=event_label),
            vjust = 0.5) +
  ggtitle("Beaver Creek 2015\n Juvenile Coho Diet (% Dry Mass)")

#### Chinook ####

#Create faceted bar charts of dry mass proportions by species and sampling event
coarselevel15 + geom_bar(aes(y=pct, x=spp, fill=TA, order=TA), 
                       data = dm_sum_final_beaver_chinook_2015,
                       stat="identity") +
  geom_text(data=dm_sum_final_beaver_chinook_2015, 
            aes(x=1, y=4, label=event_label),
            vjust = 0.5) +
  ggtitle("Beaver Creek 2015\n Juvenile Chinook Diet (% Dry Mass)")

#################
# Russian River #
#################

#### Coho ####

#Create faceted bar charts of dry mass proportions by species and sampling event
coarselevel15 + geom_bar(aes(y=pct, x=spp, fill=TA, order=TA), 
                         data = dm_sum_final_russian_coho_2015,
                         stat="identity") +
  geom_text(data=dm_sum_final_russian_coho_2015, 
            aes(x=1, y=4, label=event_label),
            vjust = 0.5) +
  ggtitle("Russian River 2015\n Juvenile Coho Diet (% Dry Mass)")

#### Chinook ####

#Create faceted bar charts of dry mass proportions by species and sampling event
coarselevel15 + geom_bar(aes(y=pct, x=spp, fill=TA, order=TA), 
                         data = dm_sum_final_russian_chinook_2015,
                         stat="identity") +
  geom_text(data=dm_sum_final_russian_chinook_2015, 
            aes(x=1, y=4, label=event_label),
            vjust = 0.5) +
  ggtitle("Russian River 2015\n Juvenile Chinook Diet (% Dry Mass)")

###################
# Ptarmigan Creek #
###################


#### Coho ####

#Create faceted bar charts of dry mass proportions by species and sampling event
coarselevel15 + geom_bar(aes(y=pct, x=spp, fill=TA, order=TA), 
                         data = dm_sum_final_ptarmigan_coho_2015,
                         stat="identity") +
  geom_text(data=dm_sum_final_ptarmigan_coho_2015, 
            aes(x=1, y=4, label=event_label),
            vjust = 0.5) +
  ggtitle("Ptarmigan Creek 2015\n Juvenile Coho Diet (% Dry Mass)")

#### Chinook ####

#Create faceted bar charts of dry mass proportions by species and sampling event
coarselevel15 + geom_bar(aes(y=pct, x=spp, fill=TA, order=TA), 
                         data = dm_sum_final_ptarmigan_chinook_2015,
                         stat="identity") +
  geom_text(data=dm_sum_final_ptarmigan_chinook_2015, 
            aes(x=1, y=4, label=event_label),
            vjust = 0.5) +
  ggtitle("Ptarmigan Creek 2015\n Juvenile Chinook Diet (% Dry Mass)")

###############
# Kenai River #
###############


#### Coho ####

#Create faceted bar charts of dry mass proportions by species and sampling event
coarselevel15 + geom_bar(aes(y=pct, x=spp, fill=TA, order=TA), 
                         data = dm_sum_final_kenai_coho_2015,
                         stat="identity") +
  geom_text(data=dm_sum_final_kenai_coho_2015, 
            aes(x=1, y=4, label=event_label),
            vjust = 0.5) +
  ggtitle("Kenai River 2015\n Juvenile Coho Diet (% Dry Mass)")

#### Chinook ####

#Create faceted bar charts of dry mass proportions by species and sampling event
coarselevel15 + geom_bar(aes(y=pct, x=spp, fill=TA, order=TA), 
                         data = dm_sum_final_kenai_chinook_2015,
                         stat="identity") +
  geom_text(data=dm_sum_final_kenai_chinook_2015, 
            aes(x=1, y=4, label=event_label),
            vjust = 0.5) +
  ggtitle("Kenai River 2015\n Juvenile Chinook Diet (% Dry Mass)")

# ------------------------------------------------------------------------ #






#######################################
#
# c.) Diet mass per fish mass by season (very coarse gut "fullness" visualization)
# 
# e.g., Terrestrial vs. Aquatic source
#
#######################################

# join T/A/M diet dry mass diet proportions of each fish to individual fish mass data.

# read in desired worksheet of diet data
diet_ta <- SCTC2015 %>% 
  gs_read(ws = "diet_T_A_M",col_names=TRUE) %>%
  
  #rename needed columns  
  rename(sample.id=X1,
         UNK=X2,
         total_mass=`Grand Total`) %>%
  
  # select needed columns
  select(sample.event,sample.id,UNK,A,M,`T`,total_mass) %>%
  
  #transform columns to desired classes
  transform(sample.id=as.factor(sample.id)) %>%
  
  #filter out NA site.id's
  filter(!is.na(sample.id),
         !is.na(sample.event),
         !is.na(total_mass)) 

# read in fish length and weight data
lw <- SCTC2015 %>% 
  gs_read(ws = "C 2015 Diet & LW Data",col_names=TRUE) %>%
  rename(wt=Weight_g,len=Length_mm) %>%
  select(river,sample.event.num,sample.event,sample.id,wt,len,spp) %>%
  filter(!is.na(sample.id),
         !is.na(sample.event),
         !is.na(wt))

# join diet data with fish length/weight data.  
fnl <- left_join(lw,diet_ta) %>%
# Mutate new column for grams of food per grams of fish.
  mutate(dm_per_fish = total_mass/wt ) %>%
  filter(total_mass != 0) %>% #considering NON-EMPTY fish only
  filter(sample.id != "LBC-F1-20150605-D11") %>% # this sample is a HUGE outlier; piscivory event
  mutate(Season = fct_recode(as.factor(sample.event.num),         #code sampling events as seasons
                             "Early Summer"="1", 
                             "Mid-Summer" = "2",
                             "Late Summer" ="3"))


# STUCK ON ERROR 8/17/17
# Spent a long time here attempting to create a theme object that could easily be 
# propograted throughout many plots for diet "fullness."  When I try to create this 
# theme object, something about the xlab, ylab, and ggtitle commands don't jive.  
# Something about applying a non-binary command to non-binary object error repeatedly
# appears despite many rearrangements.

# for now, repeat code in each iteration... argh.

# create object for each drainage
fnl_beaver <- fnl %>%
  filter(river == "Beaver Creek")
fnl_russian <- fnl %>%
  filter(river == "Russian River")
fnl_ptarmigan <- fnl %>%
  filter(river == "Ptarmigan Creek")
fnl_kenai <- fnl %>%
  filter(river == "Kenai River")

# Plot total dry mass in diet per grams of fish for each event.

# 2015 beaver creek g diet mass per g fish
bc15 <- ggplot(fnl_beaver,aes(x=factor(Season),y=dm_per_fish,color=Season)) + 
  geom_jitter(size=4) +
  theme_bw() +
  theme(plot.title = element_text(size = 26, face = "bold"),
        axis.text.x = element_text(size=14, face="bold"),
        axis.text.y = element_text(size=16, face="bold"),
        axis.title.x = element_text(size=22, face="bold"),
        axis.title.y = element_text(size = 22, face = "bold")) +
  theme(legend.title = element_text(size=30, face="bold")) + 
  theme(legend.text = element_text(size = 16, face = "bold")) + 
  ggtitle ("Beaver Creek 2015\nDiet Mass per Fish Mass ") +
  xlab(" ") +  
  ylab ("mg Diet Mass / g Fish Mass") 

# plot
bc15

# plot 2015 russian river g diet mass per g fish
rr15 <- ggplot(fnl_russian,aes(x=factor(Season),y=dm_per_fish,color=Season)) + 
  geom_jitter(size=4) +
  theme_bw() +
  theme(plot.title = element_text(size = 26, face = "bold"),
        axis.text.x = element_text(size=14, face="bold"),
        axis.text.y = element_text(size=16, face="bold"),
        axis.title.x = element_text(size=22, face="bold"),
        axis.title.y = element_text(size = 22, face = "bold")) +
  theme(legend.title = element_text(size=30, face="bold")) + 
  theme(legend.text = element_text(size = 16, face = "bold")) + 
  ggtitle ("Russian River 2015\nDiet Mass per Fish Mass ") +
  xlab(" ") +  
  ylab ("mg Diet Mass / g Fish Mass") 

# plot
rr15

# plot 2015 ptarmigan creek g diet mass per g fish
pc15 <- ggplot(fnl_ptarmigan,aes(x=factor(Season),y=dm_per_fish,color=Season)) + 
  geom_jitter(size=4) +
  theme_bw() +
  theme(plot.title = element_text(size = 26, face = "bold"),
        axis.text.x = element_text(size=14, face="bold"),
        axis.text.y = element_text(size=16, face="bold"),
        axis.title.x = element_text(size=22, face="bold"),
        axis.title.y = element_text(size = 22, face = "bold")) +
  theme(legend.title = element_text(size=30, face="bold")) + 
  theme(legend.text = element_text(size = 16, face = "bold")) + 
  ggtitle ("Ptarmigan Creek 2015\nDiet Mass per Fish Mass ") +
  xlab(" ") +  
  ylab ("mg Diet Mass / g Fish Mass") 

# plot
pc15

# plot 2015 kenai river g diet mass per g fish
kr15 <- ggplot(fnl_kenai,aes(x=factor(Season),y=dm_per_fish,color=Season)) + 
  geom_jitter(size=4) +
  theme_bw() +
  theme(plot.title = element_text(size = 26, face = "bold"),
        axis.text.x = element_text(size=14, face="bold"),
        axis.text.y = element_text(size=16, face="bold"),
        axis.title.x = element_text(size=22, face="bold"),
        axis.title.y = element_text(size = 22, face = "bold")) +
  theme(legend.title = element_text(size=30, face="bold")) + 
  theme(legend.text = element_text(size = 16, face = "bold")) + 
  ggtitle ("Kenai River 2015\nDiet Mass per Fish Mass ") +
  xlab(" ") +  
  ylab ("mg Diet Mass / g Fish Mass") 

# plot
kr15