###########################################################################
###########################################################################
###The Uncertain Bioenergetics of North Atlantic right whales (R Script)################
###########################################################################
###########################################################################
#Authors: Jasmin C. HŸtt, Peter Corkeron, Julie M. van der Hoop, and Michael J. Moore
#Corresponding author email: JasHuett@gmail.com
#Script published as Electronic Supplementary Material

###Basic R Settings###
rm (list = ls (all = TRUE)) #Clear Workspace
setwd("D:/") #Insert your working directory here
set.seed(12345)
#load libraries
require('ggplot2')
library(ggpubr)
library(ggprism)
require('tidyverse')
require('ggdist') #for raincloud plots
#Set up ggplot theme
theme_custom <- function (base_size = 11, base_family = "") {
theme_classic() %+replace%
theme(
axis.line = element_line(),
axis.ticks = element_line(),
axis.text = element_text(),
text = element_text(size = 9),
legend.position="none"
)
}
theme_set(theme_custom())
#######################################################
###First phase of this work - food and energy income###
#######################################################
###Parameters###
#Morphological and general ones
L<-c(1293.67,1327.0,1362.57) #Body Length (cm) from Christiansen et al 2020 and Miller et al 2020 for NARWnow, NARW2000, and
SRW
HW<-c(180.0298,202.0,179.3783) #Head Width (10% BL from rostrum, cm) from Christiansen et al 2020 and Miller et al 2020 for
NARWnow, NARW2000, and SRW
BAL<-(0.2077*(L/100)-1.095) ##Baleen Length (m) for NARWnow, NARW2000, and SRW
G<-((HW/100)*BAL)/2 #Gape Area (m^2) for NARWnow, NARW2000, and SRW
L<-append(L,1050) #Add Body Length (cm) for whale 3617 from Stewart et al 2021
G<-append(G,0.9) #Add Gape Area (m^2) for whale 3617 with body length from Stewart et al 2021 but Gape read from Figure 2 in
Van der Hoop et al 2019
Vol<-c(23.21,30.04,34.45,14.845) #Body Volume (m^3) from Christiansen et al 2020, Miller et al 2020, and Stewart et al 2021
M<-754.63*Vol ##Mass (kg)
#Filtration estimate
S<-(0.09*L)/100 #Swimming Speed (m/s)
Slurped<-S*G #Amount of water swept by a whale
EFc<-runif(10000, min=0.885, max=1.5) #Range of capture efficiency
#Calculate seawater filtration rate (m^3/s)
#Use outer to combine Slurped & EFc for all cases of both
Filtered<-outer(EFc,Slurped,FUN='*') #Water swept by a whale (m^3/s)
#########################################################
#Foraging time
#Fa in hours ranges from 15.8 to 22
#Time spent foraging
ForageHrs<-runif(10000, min=15.8,max=22) #From Kenney 1986 and Fortune et al. 2013
Fa<-ForageHrs*60*60 #Time spent for foraging activities (s/day)
###########################################################
#Propoortion of the dives spent feeding from Baumgartner & Mate 2003, Table 3
MarksTime<-as.data.frame(read.table('MarksTimeFeeding.txt'))
names(MarksTime)<-'Proportion'
str(MarksTime)
hist(MarksTime$Proportion)
#With Log-Normal distribution sampling of Calanus abundance
#Code to produce the log-normal distribution lifted from:
#https://gist.github.com/msalganik/7d8762e69b3954b3a5a2ec437f62d6e9
#and modded for here
m<-mean(MarksTime$Proportion)/100
s<-sd(MarksTime$Proportion)/100
location <- log(m^2 / sqrt(s^2 + m^2))
shape <- sqrt(log(1 + (s^2 / m^2)))
print(paste("location:", location))
print(paste("shape:", shape))
#Checking this works
draws3 <- rlnorm(n=100000,location,shape)
mean(draws3)
sd(draws3) #it works!
#Log-Normal distribution
Fe<-rlnorm(n=10000,location, shape)
Ft<-Fa*Fe ##Time spent feeding effectively (s/day)
#########################################################
#Water filtered per day
DailyWF<-as.data.frame(Filtered*Ft) #m^3 water filtered/ind./day
names(DailyWF)<-c('NARWnow','NARW2000','SRW',"Stunted")
DailyWF2<-stack(DailyWF)
#Divide values by 1000 to be more readable
DailyWF2$values<-DailyWF2$values/1000
#Raincloud plots of each of these
#Code for Raincloud plots taken from https://kirikuroda.com/en/post/2021/10/22/raincloudplot-ggplot2/
p1_r <- ggplot(DailyWF2, aes(x=values, y=ind, color=ind,
fill=ind)) +
stat_halfeye(
point_color=NA, .width=0, height=0.6,
position=position_nudge(y = 0.2)
) +
geom_boxplot(
position=position_nudge(y = 0.1),
width=0.1, outlier.shape=NA, alpha=0.1
)+
xlim(0,220)+
labs(x=expression("Water filtered [1000 m"^"3"*" day"^"-1"*"]"))+
scale_colour_manual(values = c("black", "black","black","black"))+
scale_fill_manual(values = c("black", "black","black","black"))+
theme(axis.title.y= element_blank())
p1_r
#Values for only the three classes
DailyWF2Classes <- DailyWF2[DailyWF2$ind=="NARWnow"|DailyWF2$ind=="NARW2000"|DailyWF2$ind=="SRW",]
mean(DailyWF2Classes$values)*1000
sd(DailyWF2Classes$values)*1000
#95% bounds
#sort in order
DailyWF2Classes.sorted<-DailyWF2Classes[order(DailyWF2Classes$values),]
#5%
round(DailyWF2Classes.sorted[1500,1],3)
#95%
round(DailyWF2Classes.sorted[28500,1],3)
#Values all
mean(DailyWF2$values)*1000
sd(DailyWF2$values)*1000
#95% bounds
#sort in order
DailyWF2.sorted<-DailyWF2[order(DailyWF2$values),]
#5%
round(DailyWF2.sorted[2000,1],3)
#95%
round(DailyWF2.sorted[40000,1],3)
#by type or right whale
round(colMeans(DailyWF),0)
round(apply(DailyWF,2,sd),1)
########################################################
#How much plankton?
CFmass<-0.25 #Dry mass of C. finmarchicus (mg/individual) - average from Maps et al 2010, and Maps et al 2011
CFmassw<-9.44061*CFmass^0.946 #Wet mass of C. finmarchicus (mg/individual) - average - conversion after Wiebe et al 1988
#Plankton density (ind./m^3)from Baumgartner & Mate MEPS 2003 Table 4
MarksData<-as.data.frame(read.table('MarksCopepods.txt'))
#get it in order
MarksData<-as.data.frame(sort(MarksData$V1))
names(MarksData)<-'Density'
mean(MarksData$Density) #just checking!
sd(MarksData$Density)
#Let's look at their data
ggplot(MarksData, aes(x=Density))+geom_dotplot()
#There's not that much there.
###################################################
#Resample plankton data
#With Log-Normal distribution sampling of Calanus abundance
#Code to produce the log-normal distribution lifted from:
#https://gist.github.com/msalganik/7d8762e69b3954b3a5a2ec437f62d6e9
#and modded for here
m<-mean(MarksData$Density)
s<-sd(MarksData$Density)
location <- log(m^2 / sqrt(s^2 + m^2))
shape <- sqrt(log(1 + (s^2 / m^2)))
print(paste("location:", location))
print(paste("shape:", shape))
#Checking this works
draws3 <- rlnorm(n=100000,location,shape)
mean(draws3)
sd(draws3) #it works!
#Log-Normal distribution
Cd<-rlnorm(n=10000,location, shape)
#Income with average food and dry weight
FDmg<-Filtered*Cd*CFmass*Ft ##Food Income (mg/day)
#Being clear on the division from mg to tonnes
FDkg<-FDmg/(10^6) ##Food Income (kg/day)
FDtonne<-FDkg/1000 ##Food Income (tonnes/day)
##see how the kg/day data look for 6K plankton/m^3
FDkg<-as.data.frame(FDkg)
names(FDkg)<-c('NARWnow','NARW2000','SRW',"Stunted")
FDkg2<-stack(FDkg)
#see how this looks
p3_r <- ggplot(FDkg2, aes(x=values, y=ind, color=ind,
fill=ind)) +
stat_halfeye(
point_color=NA, .width=0, height=0.6,
position=position_nudge(y = 0.2)
) +
geom_boxplot(
position=position_nudge(y = 0.1),
width=0.1, outlier.shape=NA, alpha=0.1
)+
xlim(0,500)+
labs(x=expression("Food income [kg d"^"-1"*"]"))+
theme(axis.title.y= element_blank())+
scale_colour_manual(values = c("black", "black","black","black"))+
scale_fill_manual(values = c("black", "black","black","black"))
p3_r
#########################
#Means and SDs for only the three classes
FDkg2Classes <- FDkg2[FDkg2$ind=="NARW2000"|FDkg2$ind=="NARWnow"|FDkg2$ind=="SRW",]
round(mean(FDkg2Classes$values),0)
round(sd(FDkg2Classes$values),1)
#95% bounds
#sort in order
FDkg2Classes.sorted<-FDkg2Classes[order(FDkg2Classes$values),]
#5%
round(FDkg2Classes.sorted[1500,1],3)
#95%
round(FDkg2Classes.sorted[28500,1],3)
#Estimates all
round(mean(FDkg2$values),0)
round(sd(FDkg2$values),1)
#95% bounds
#sort in order
FDkg2.sorted<-FDkg2[order(FDkg2$values),]
#5%
round(FDkg2.sorted[2000,1],3)
#95%
round(FDkg2.sorted[38000,1],3)
#
round(min(FDkg2$values),1)
round(max(FDkg2$values),1)
round(colMeans(FDkg),0)
round(apply(FDkg,2,sd),1)
#################################################
#Energy in each plankter - estimates
#Michaud & Taggart 2007: 3.4 - 5J/ind
#McKinstry et al 2013 Table 2: 8.02-12.67 (C5ec)
Ev<-runif(10000, min=3.4, max=12.67) #Range of energy content per plankter
#Income with average food and plankton energy content
Jin<-Filtered*Ev*Cd*Ft
##see how the J/day data look for 6K plankton/m^3
Jin<-as.data.frame(Jin)
names(Jin)<-c('NARWnow','NARW2000','SRW',"Stunted")
Jin2<-stack(Jin)
Jin2$MJ<-Jin2$values/(10^6)
#see how this looks
p4_r <- ggplot(Jin2, aes(x=MJ, y=ind, color=ind,
fill=ind)) +
stat_halfeye(
point_color=NA, .width=0, height=0.6,
position=position_nudge(y = 0.2)
) +
geom_boxplot(
position=position_nudge(y = 0.1),
width=0.1, outlier.shape=NA, alpha=0.1
)+
xlim(0,20000)+
labs(x=expression("Energy income [MJ d"^"-1"*"]"))+
theme(axis.title.y= element_blank())+
scale_colour_manual(values = c("black", "black","black","black"))+
scale_fill_manual(values = c("black", "black","black","black"))
p4_r
#Estimates for only the three classes
Jin2Classes <- Jin2[Jin2$ind=="NARWnow"|Jin2$ind=="NARW2000"|Jin2$ind=="SRW",]
round(mean(Jin2Classes$MJ),0)
round(sd(Jin2Classes$MJ),1)
#95% bounds
#sort in order
Jin2Classes.sorted<-Jin2Classes[order(Jin2Classes$MJ),]
#5%
round(Jin2Classes.sorted[1500,1]/10^6,2)
#95%
round(Jin2Classes.sorted[28500,1]/10^6,2)
#Estimates all
round(mean(Jin2$MJ),0)
round(sd(Jin2$MJ),1)
#95% bounds
#sort in order
Jin2.sorted<-Jin2[order(Jin2$MJ),]
#5%
round(Jin2.sorted[2000,1]/10^6,2)
#95%
round(Jin2.sorted[38000,1]/10^6,2)
#
round(min(Jin2$MJ),1)
round(max(Jin2$MJ),1)
round(colMeans(Jin)/10^6,0)
round(apply(Jin,2,sd)/10^6,1)
####################################################
###Second phase of this work - energy expenditure###
####################################################
#BMR calculations
#Kleiber's simple value
KBMR<-292.88*(M^0.75) #Kleiber estimate in kilojoules/day
KBMR
KBMR.Mj<-KBMR/1000
round(KBMR.Mj,0) #Kleiber estimate in megaJoules/day
#for ease of viewing
######################################################
#BMR range as 4 values
BMRrange<-c(0.3,0.5,0.75,1.0)
#range of potential values for daily BMR as proportion of Kleiber
BMRdaily<-outer(BMRrange,KBMR.Mj,FUN='*') #range of estimates of daily BMR
#What does this look like?
BMRdaily<-as.data.frame(BMRdaily)
names(BMRdaily)<-c('NARWnow','NARW2000','SRW',"Stunted")
BMRdaily2<-stack(BMRdaily)
BMRdaily2<-cbind(BMRdaily2,BMRrange)
names(BMRdaily2)<-c('BMRinMJ','Class','BMRMultiplier')
BMRdaily2$Class<-as.factor(BMRdaily2$Class)
BMRdaily2$BMRMultiplier<-as.factor(BMRdaily2$BMRMultiplier)
#see how this looks
#Very basic bar graph
p5<-ggplot(data=BMRdaily2, aes(x=BMRinMJ, y=Class,fill=BMRMultiplier)) +
geom_bar(stat="identity", colour="black", position=position_dodge())+
theme(axis.title.y= element_blank())+
labs(x=expression("BMR [MJ d"^"-1"*"]"))+
theme(legend.position = "right")+
guides(fill=guide_legend(title="Factor of Kleiber"))
p5
#######################################
#Estimates for only the three classes
round(mean(BMRdaily2$BMRinMJ[BMRdaily2$Class=="NARWnow"|BMRdaily2$Class=="NARW2000"|BMRdaily2$Class=="SRW"]),1)
round(sd(BMRdaily2$BMRinMJ[BMRdaily2$Class=="NARWnow"|BMRdaily2$Class=="NARW2000"|BMRdaily2$Class=="SRW"]),2)
round(min(BMRdaily2$BMRinMJ[BMRdaily2$Class=="NARWnow"|BMRdaily2$Class=="NARW2000"|BMRdaily2$Class=="SRW"]),1)
round(max(BMRdaily2$BMRinMJ[BMRdaily2$Class=="NARWnow"|BMRdaily2$Class=="NARW2000"|BMRdaily2$Class=="SRW"]),1)
#Estimates for specific factors of Kleiber
BMRdaily2Classes <- BMRdaily2[BMRdaily2$Class=="NARWnow"|BMRdaily2$Class=="NARW2000"|BMRdaily2$Class=="SRW",]
mean(BMRdaily2Classes$BMRinMJ[BMRdaily2Classes$BMRMultiplier=="0.3"])
sd(BMRdaily2Classes$BMRinMJ[BMRdaily2Classes$BMRMultiplier=="0.3"])
mean(BMRdaily2Classes$BMRinMJ[BMRdaily2Classes$BMRMultiplier=="1"])
sd(BMRdaily2Classes$BMRinMJ[BMRdaily2Classes$BMRMultiplier=="1"])
#No CIs here as there are only 3 values
###############################################
#Compare daily BMR with daily energy income###
###############################################
#Daily BMR with daily energy income and allowing for BMR sensitivity analysis
#sapply used here
NARWnow<-as.data.frame(t(round(sapply(Jin[,1], function(i)i/BMRdaily[,1]/10^6),2)))
names(NARWnow)<-c('BMR*0.3','BMR*0.5','BMR*0.75','BMR*1')
NARW2000<-as.data.frame(t(round(sapply(Jin[,2], function(i)i/BMRdaily[,2]/10^6),2)))
names(NARW2000)<-c('BMR*0.3','BMR*0.5','BMR*0.75','BMR*1')
SRW<-as.data.frame(t(round(sapply(Jin[,3], function(i)i/BMRdaily[,3]/10^6),2)))
names(SRW)<-c('BMR*0.3','BMR*0.5','BMR*0.75','BMR*1')
Stunted<-as.data.frame(t(round(sapply(Jin[,4], function(i)i/BMRdaily[,4]/10^6),2)))
names(Stunted)<-c('BMR*0.3','BMR*0.5','BMR*0.75','BMR*1')
NARWnow.2<-stack(NARWnow)
NARW2000.2<-stack(NARW2000)
SRW.2<-stack(SRW)
Stunted.2<-stack(Stunted)
NARWnow.2<-cbind(NARWnow.2, 'NARWnow')
NARW2000.2<-cbind(NARW2000.2,'NARW2000')
SRW.2<-cbind(SRW.2, 'SRW')
Stunted.2<-cbind(Stunted.2, 'Stunted')
names(NARWnow.2)<-c('Income','BMRMultiplier','Class')
names(NARW2000.2)<-c('Income','BMRMultiplier','Class')
names(SRW.2)<-c('Income','BMRMultiplier','Class')
names(Stunted.2)<-c('Income','BMRMultiplier','Class')
#############################################
p7<-ggplot(NARWnow.2, aes(x=Income, y=BMRMultiplier, color=BMRMultiplier,
fill=BMRMultiplier)) +
stat_halfeye(
point_color=NA, .width=0, height=0.6,
position=position_nudge(y = 0.3)
) +
geom_boxplot(
position=position_nudge(y = 0.2),
width=0.1, outlier.shape=NA, alpha=0.1
)+
xlim(0,80)+
labs(x=expression("Energy income as a multiplier of BMR: NARW now"))
p7
p8<-ggplot(NARW2000.2, aes(x=Income, y=BMRMultiplier, color=BMRMultiplier,
fill=BMRMultiplier)) +
stat_halfeye(
point_color=NA, .width=0, height=0.6,
position=position_nudge(y = 0.3)
) +
geom_boxplot(
position=position_nudge(y = 0.2),
width=0.1, outlier.shape=NA, alpha=0.1
)+
xlim(0,80)+
labs(x=expression("Energy income as a multiplier of BMR: NARW2000"))
p8
p9<-ggplot(SRW.2, aes(x=Income, y=BMRMultiplier, color=BMRMultiplier,
fill=BMRMultiplier)) +
stat_halfeye(
point_color=NA, .width=0, height=0.6,
position=position_nudge(y = 0.3)
) +
geom_boxplot(
position=position_nudge(y = 0.2),
width=0.1, outlier.shape=NA, alpha=0.1
)+
xlim(0,80)+
labs(x=expression("Energy income as a multiplier of BMR: SRW"))
p9
ps<-ggplot(Stunted.2, aes(x=Income, y=BMRMultiplier, color=BMRMultiplier,
fill=BMRMultiplier)) +
stat_halfeye(
point_color=NA, .width=0, height=0.6,
position=position_nudge(y = 0.3)
) +
geom_boxplot(
position=position_nudge(y = 0.2),
width=0.1, outlier.shape=NA, alpha=0.1
)+
xlim(0,80)+
labs(x=expression("Energy income as a multiplier of BMR: Stunted"))
ps
#Plot all together
#join them all up
RW.all<-rbind(NARWnow.2,NARW2000.2,SRW.2,Stunted.2)
#Create a final category
RW.all$Class2<-paste(RW.all$BMRMultiplier,RW.all$Class, sep="")
p10<-ggplot(RW.all, aes(x=Income, y=Class2, color=Class2,
fill=Class2)) +
stat_halfeye(
point_color=NA, .width=0, height=0.6,
position=position_nudge(y = 0.3)
) +
geom_boxplot(
position=position_nudge(y = 0.2),
width=0.1, outlier.shape=NA, alpha=0.1
)+
xlim(0,80)+
labs(x=expression("Energy income as a multiplier of BMR"))+
theme(axis.title.y= element_blank())+
scale_colour_manual(values = c("#F8766D","#F8766D","#F8766D","#F8766D","#7CAE00","#7CAE00","#7CAE00","#7CAE00"
,"#00BFC4","#00BFC4","#00BFC4","#00BFC4","#C77CFF","#C77CFF","#C77CFF","#C77CFF"))+
scale_fill_manual(values = c("#F8766D","#F8766D","#F8766D","#F8766D","#7CAE00","#7CAE00","#7CAE00","#7CAE00"
,"#00BFC4","#00BFC4","#00BFC4","#00BFC4","#C77CFF","#C77CFF","#C77CFF","#C77CFF"))+
scale_y_discrete(limits=rev,labels =
c("Stunted","SRW","NARWnow","NARW2000","Stunted","SRW","NARWnow","NARW2000","Stunted","SRW","NARWnow","NARW2000","
Stunted","SRW","NARWnow","NARW2000"))
p10
######################################################
#What do these values look like?
tapply(RW.all$Income, RW.all$Class2, summary)
RW.all.classes <- RW.all[RW.all$Class=="NARWnow"|RW.all$Class=="NARW2000"|RW.all$Class=="SRW",]
mean(RW.all.classes$Income[RW.all.classes$BMRMultiplier=="BMR*0.3"])
sd(RW.all.classes$Income[RW.all.classes$BMRMultiplier=="BMR*0.3"])
mean(RW.all.classes$Income[RW.all.classes$BMRMultiplier=="BMR*1"])
sd(RW.all.classes$Income[RW.all.classes$BMRMultiplier=="BMR*1"])
#No CIs here as there are only three values included
############################################
#Energetic costs for drag of entanglement###
############################################
#from van der Hoop et al 2017
#range 7.24*10^7 to 7.52*10^8 J/day
#convert to MJ/day for comparison: divide by 10^6
#########################
names(Jin2)<-c("kJintake","Class","MJintake")
Tanglerange<-runif(10000, min=72.4,max=752) #in MJ
###############################################
#Daily entanglement cost relative to food intake
Jin.Tngl<-cbind(Jin2,Tanglerange) #in MJ
Jin.Tngl$TangleToIntake<-Jin.Tngl$Tanglerange/Jin.Tngl$MJintake
#Calculate Entanglement costs
mean(Jin.Tngl$Tanglerange)
sd(Jin.Tngl$Tanglerange)
#plot this
p11<-ggplot(Jin.Tngl, aes(x=TangleToIntake, y=Class, color=Class,
fill=Class)) +
stat_halfeye(
point_color=NA, .width=0, height=0.6,
position=position_nudge(y = 0.2)
) +
geom_boxplot(
position=position_nudge(y = 0.1),
width=0.1, outlier.shape=NA, alpha=0.1
)+
xlim(0,1.5)+
labs(x=expression("Entanglement expenditure per energy income"))+
scale_colour_manual(values = c("black", "black","black","black"))+
scale_fill_manual(values = c("black", "black","black","black"))+
theme(axis.title.y= element_blank())
p11
############################################
#Estimates for only the three classes
Jin.TnglClasses <- Jin.Tngl[Jin.Tngl$Class=="NARWnow"|Jin.Tngl$Class=="NARW2000"|Jin.Tngl$Class=="SRW",]
round(mean(Jin.TnglClasses$TangleToIntake),3)*100
round(sd(Jin.TnglClasses$TangleToIntake),4)*100
#
#95% CIs
Jin.TnglClasses.sorted<-Jin.TnglClasses[order(Jin.TnglClasses$TangleToIntake),]
#5%
round(Jin.TnglClasses.sorted[1500,5],4)*100
#95%
round(Jin.TnglClasses.sorted[28500,5],4)*100
#These estimates as percentages for all
round(mean(Jin.Tngl$TangleToIntake),3)*100
round(sd(Jin.Tngl$TangleToIntake),4)*100
#
#95% CIs
Jin.Tngl.sorted<-Jin.Tngl[order(Jin.Tngl$TangleToIntake),]
#5%
round(Jin.Tngl.sorted[2000,5],4)*100
#95%
round(Jin.Tngl.sorted[38000,5],4)*100
#
round(min(Jin.Tngl$TangleToIntake),3)*100
round(max(Jin.Tngl$TangleToIntake),3)*100
aggregate(Jin.Tngl$TangleToIntake, list(Jin.Tngl$Class), FUN=mean)
aggregate(Jin.Tngl$TangleToIntake, list(Jin.Tngl$Class), FUN=sd)
############################################
#Daily entanglement cost relative to BMR###
###########################################
NARWnow.Tngl<-as.data.frame(t(round(sapply(Tanglerange, function(i)i/BMRdaily[,1]),3)))
names(NARWnow.Tngl)<-c('BMR*0.3','BMR*0.5','BMR*0.75','BMR*1')
NARW2000.Tngl<-as.data.frame(t(round(sapply(Tanglerange, function(i)i/BMRdaily[,2]),3)))
names(NARW2000.Tngl)<-c('BMR*0.3','BMR*0.5','BMR*0.75','BMR*1')
SRW.Tngl<-as.data.frame(t(round(sapply(Tanglerange, function(i)i/BMRdaily[,3]),3)))
names(SRW.Tngl)<-c('BMR*0.3','BMR*0.5','BMR*0.75','BMR*1')
Stunted.Tngl<-as.data.frame(t(round(sapply(Tanglerange, function(i)i/BMRdaily[,4]),3)))
names(Stunted.Tngl)<-c('BMR*0.3','BMR*0.5','BMR*0.75','BMR*1')
NARWnow.Tngl.2<-stack(NARWnow.Tngl)
NARW2000.Tngl.2<-stack(NARW2000.Tngl)
SRW.Tngl.2<-stack(SRW.Tngl)
Stunted.Tngl.2<-stack(Stunted.Tngl)
NARWnow.Tngl.2<-cbind(NARWnow.Tngl.2, 'NARWnow')
NARW2000.Tngl.2<-cbind(NARW2000.Tngl.2,'NARW2000')
SRW.Tngl.2<-cbind(SRW.Tngl.2, 'SRW')
Stunted.Tngl.2<-cbind(Stunted.Tngl.2, 'Stunted')
names(NARWnow.Tngl.2)<-c('TangleCost','BMRMultiplier','Class')
names(NARW2000.Tngl.2)<-c('TangleCost','BMRMultiplier','Class')
names(SRW.Tngl.2)<-c('TangleCost','BMRMultiplier','Class')
names(Stunted.Tngl.2)<-c('TangleCost','BMRMultiplier','Class')
#################
#Join all of these
RW.Tngl.all<-rbind(NARWnow.Tngl.2,NARW2000.Tngl.2,SRW.Tngl.2, Stunted.Tngl.2)
#Create a final category
RW.Tngl.all$Class2<-paste(RW.all$BMRMultiplier,RW.all$Class, sep="")
#plot this
p12<-ggplot(RW.Tngl.all, aes(x=TangleCost, y=Class2, color=Class2,
fill=Class2)) +
stat_halfeye(
point_color=NA, .width=0, height=0.6,
position=position_nudge(y = 0.3)
) +
geom_boxplot(
position=position_nudge(y = 0.2),
width=0.1, outlier.shape=NA, alpha=0.1
)+
xlim(0,6)+
labs(x=expression("Entanglement expenditure per BMR"))+
theme(axis.title.y= element_blank())+
scale_colour_manual(values = c("#F8766D","#F8766D","#F8766D","#F8766D","#7CAE00","#7CAE00","#7CAE00","#7CAE00"
,"#00BFC4","#00BFC4","#00BFC4","#00BFC4","#C77CFF","#C77CFF","#C77CFF","#C77CFF"))+
scale_fill_manual(values = c("#F8766D","#F8766D","#F8766D","#F8766D","#7CAE00","#7CAE00","#7CAE00","#7CAE00"
,"#00BFC4","#00BFC4","#00BFC4","#00BFC4","#C77CFF","#C77CFF","#C77CFF","#C77CFF"))+
scale_y_discrete(limits=rev,labels =
c("Stunted","SRW","NARWnow","NARW2000","Stunted","SRW","NARWnow","NARW2000","Stunted","SRW","NARWnow","NARW2000","
Stunted","SRW","NARWnow","NARW2000"))
p12
############################################
#What do these values look like?
tapply(RW.Tngl.all$TangleCost, RW.Tngl.all$Class2, summary)
RW.Tngl.all.classes <- RW.Tngl.all[RW.Tngl.all$Class=="NARWnow"|RW.Tngl.all$Class=="NARW2000"|RW.Tngl.all$Class=="SRW",]
mean(RW.Tngl.all.classes$TangleCost[RW.Tngl.all.classes$BMRMultiplier=="BMR*0.3"])
sd(RW.Tngl.all.classes$TangleCost[RW.Tngl.all.classes$BMRMultiplier=="BMR*0.3"])
mean(RW.Tngl.all.classes$TangleCost[RW.Tngl.all.classes$BMRMultiplier=="BMR*1"])
sd(RW.Tngl.all.classes$TangleCost[RW.Tngl.all.classes$BMRMultiplier=="BMR*1"])
mean(RW.Tngl.all.classes$TangleCost[RW.Tngl.all.classes$BMRMultiplier=="BMR*0.5"])
sd(RW.Tngl.all.classes$TangleCost[RW.Tngl.all.classes$BMRMultiplier=="BMR*0.5"])
###########################################################################
###########################################################################
###Now let's assume - third phase - that entanglement affects filtration###
###########################################################################
###########################################################################
#Filtration estimate
EFc.ent<-runif(10000, min=0.4425, max=0.75) #Range of capture efficiency
#Assuming efficiency can be halved by entanglement
#calculate seawater filtration rate (m^3/s)
#Use outer to combine Slurped & EFc for all cases of both
Filtered.ent<-outer(EFc.ent,Slurped,FUN='*') #water swept by a whale (m^3/s)
#########################################################
#Water filtered per day
DailyWF.ent<-as.data.frame(Filtered.ent*Ft) #m^3 water filtered/ind./day
names(DailyWF.ent)<-c('NARWnow','NARW2000','SRW', "Stunted")
DailyWF.ent2<-stack(DailyWF.ent)
#Divide values by 1000 to be more readable
DailyWF.ent2$values<-DailyWF.ent2$values/1000
#plot
p13_r<-ggplot(DailyWF.ent2, aes(x=values, y=ind, color=ind,
fill=ind)) +
stat_halfeye(
point_color=NA, .width=0, height=0.6,
position=position_nudge(y = 0.3)
) +
geom_boxplot(
position=position_nudge(y = 0.2),
width=0.1, outlier.shape=NA, alpha=0.1
)+
xlim(0,130)+
labs(x=expression("Water filtered (E) [1000 m"^"3"*" day"^"-1"*"]"))+
theme(axis.title.y= element_blank())+
scale_colour_manual(values = c("black", "black","black","black"))+
scale_fill_manual(values = c("black", "black","black","black"))+
theme(axis.title.y= element_blank(), axis.text.y=element_blank(), axis.ticks.y=element_blank())
p13_r
#Means and SDs for only the three classes
DailyWF.ent2Classes <- DailyWF.ent2[DailyWF.ent2$ind=="NARWnow"|DailyWF.ent2$ind=="NARW2000"|DailyWF.ent2$ind=="SRW",]
round(mean(DailyWF.ent2Classes$values),3)
round(sd(DailyWF.ent2Classes$values),4)
#95% CIs
DailyWF.ent2.Classessorted<-DailyWF.ent2Classes[order(DailyWF.ent2Classes$values),]
#5%
round(DailyWF.ent2.Classessorted[1500,1],3)
#95%
round(DailyWF.ent2.Classessorted[28500,1],3)
#values for all
round(mean(DailyWF.ent2$values),3)
round(sd(DailyWF.ent2$values),4)
#95% CIs
DailyWF.ent2.sorted<-DailyWF.ent2[order(DailyWF.ent2$values),]
#5%
round(DailyWF.ent2.sorted[2000,1],3)
#95%
round(DailyWF.ent2.sorted[38000,1],3)
#
round(min(DailyWF.ent2$values),3)
round(max(DailyWF.ent2$values),3)
#by type of right whale
round(colMeans(DailyWF.ent/1000),3)
round(apply(DailyWF.ent/1000,2,sd),3)
#compare to without entanglement
round(colMeans(DailyWF/1000),0)
round(apply(DailyWF/1000,2,sd),1)
##################################################
#Income with average food and dry weight
FDmg.ent<-Filtered.ent*Cd*CFmass*Ft ##Food Income (mg/day)
#Being clear on the division from mg to tonnes
FDkg.ent<-FDmg.ent/(10^6) ##Food Income (kg/day)
FDtonne.ent<-FDkg.ent/1000 ##Food Income (tonnes/day)
##see how the kg/day data look for 6K plankton/m^3
FDkg.ent<-as.data.frame(FDkg.ent)
names(FDkg.ent)<-c('NARWnow','NARW2000','SRW',"Stunted")
FDkg.ent2<-stack(FDkg.ent)
#see how this looks
p14_r<-ggplot(FDkg.ent2, aes(x=values, y=ind, color=ind,
fill=ind)) +
stat_halfeye(
point_color=NA, .width=0, height=0.6,
position=position_nudge(y = 0.3)
) +
geom_boxplot(
position=position_nudge(y = 0.2),
width=0.1, outlier.shape=NA, alpha=0.1
)+
xlim(0,300)+
labs(x=expression("Food income (E) [kg d"^"-1"*"]"))+
theme(axis.title.y= element_blank())+
theme(axis.title.y= element_blank(), axis.text.y=element_blank(), axis.ticks.y=element_blank())+
scale_colour_manual(values = c("black", "black","black","black"))+
scale_fill_manual(values = c("black", "black","black","black"))
p14_r
#########################
#Estimates for only the three classes
FDkg.ent2Classes <- FDkg.ent2[FDkg.ent2$ind=="NARWnow"|FDkg.ent2$ind=="NARW2000"|FDkg.ent2$ind=="SRW",]
round(mean(FDkg.ent2Classes$values),0)
round(sd(FDkg.ent2Classes$values),1)
#95% CIs
FDkg.ent2.sortedClasses<-FDkg.ent2Classes[order(FDkg.ent2Classes$values),]
#5%
round(FDkg.ent2.sortedClasses[1500,1],2)
#95%
round(FDkg.ent2.sortedClasses[28500,1],2)
#Estimates for all
round(mean(FDkg.ent2$values),0)
round(sd(FDkg.ent2$values),1)
#95% CIs
FDkg.ent2.sorted<-FDkg.ent2[order(FDkg.ent2$values),]
#5%
round(FDkg.ent2.sorted[2000,1],2)
#95%
round(FDkg.ent2.sorted[38000,1],2)
#
round(colMeans(FDkg.ent),0)
round(apply(FDkg.ent,2,sd),1)
#################################################
#Income with average food and plankton energy content
Jin.ent<-Filtered.ent*Ev*Cd*Ft
##see how the J/day data look for 6K plankton/m^3
Jin.ent<-as.data.frame(Jin.ent)
names(Jin.ent)<-c('NARWnow','NARW2000','SRW',"Stunted")
Jin.ent2<-stack(Jin.ent)
Jin.ent2$MJ<-Jin.ent2$values/(10^6)
#see how this looks
p15_r<-ggplot(Jin.ent2, aes(x=MJ, y=ind, color=ind,
fill=ind)) +
stat_halfeye(
point_color=NA, .width=0, height=0.6,
position=position_nudge(y = 0.3)
) +
geom_boxplot(
position=position_nudge(y = 0.2),
width=0.1, outlier.shape=NA, alpha=0.1
)+
xlim(0,15000)+
labs(x=expression("Energy income (E) [MJ d"^"-1"*"]"))+
theme(axis.title.y= element_blank())+
theme(axis.title.y= element_blank(), axis.text.y=element_blank(), axis.ticks.y=element_blank())+
scale_colour_manual(values = c("black", "black","black","black"))+
scale_fill_manual(values = c("black", "black","black","black"))
p15_r
#Estimates for only the three classes
Jin.ent2Classes <- Jin.ent2[Jin.ent2$Class=="NARWnow"|Jin.ent2$Class=="NARW2000"|Jin.ent2$Class=="SRW",]
round(mean(Jin.ent2Classes$values),0)
round(sd(Jin.ent2Classes$values),1)
#95% CIs
Jin.ent2.sortedClasses<-Jin.ent2Classes[order(Jin.ent2Classes$values),]
#5%
round(Jin.ent2.sortedClasses[1500,3],1)
#95%
round(Jin.ent2.sortedClasses[28500,3],1)
#Estimates for all
round(mean(Jin.ent2$MJ),0)
round(sd(Jin.ent2$MJ),1)
#95% CIs
Jin.ent2.sorted<-Jin.ent2[order(Jin.ent2$values),]
#5%
round(Jin.ent2.sorted[2000,3],1)
#95%
round(Jin.ent2.sorted[38000,3],1)
#
round(colMeans(Jin.ent)/10^6,0)
round(apply(Jin.ent,2,sd)/10^6,1)
############################################################
#sapply used here
NARWnow.ent<-as.data.frame(t(round(sapply(Jin.ent[,1], function(i)i/BMRdaily[,1]/10^6),2)))
names(NARWnow.ent)<-c('BMR*0.3','BMR*0.5','BMR*0.75','BMR*1')
NARW2000.ent<-as.data.frame(t(round(sapply(Jin.ent[,2], function(i)i/BMRdaily[,2]/10^6),2)))
names(NARW2000.ent)<-c('BMR*0.3','BMR*0.5','BMR*0.75','BMR*1')
SRW.ent<-as.data.frame(t(round(sapply(Jin.ent[,3], function(i)i/BMRdaily[,3]/10^6),2)))
names(SRW.ent)<-c('BMR*0.3','BMR*0.5','BMR*0.75','BMR*1')
Stunted.ent<-as.data.frame(t(round(sapply(Jin.ent[,4], function(i)i/BMRdaily[,4]/10^6),2)))
names(Stunted.ent)<-c('BMR*0.3','BMR*0.5','BMR*0.75','BMR*1')
NARWnow.ent.2<-stack(NARWnow.ent)
NARW2000.ent.2<-stack(NARW2000.ent)
SRW.ent.2<-stack(SRW.ent)
Stunted.ent.2<-stack(Stunted.ent)
NARWnow.ent.2<-cbind(NARWnow.ent.2, 'NARWnow.ent')
NARW2000.ent.2<-cbind(NARW2000.ent.2,'NARW2000.ent')
SRW.ent.2<-cbind(SRW.ent.2, 'SRW.ent')
Stunted.ent.2<-cbind(Stunted.ent.2, 'Stunted.ent')
names(NARWnow.ent.2)<-c('Income','BMRMultiplier','Class')
names(NARW2000.ent.2)<-c('Income','BMRMultiplier','Class')
names(SRW.ent.2)<-c('Income','BMRMultiplier','Class')
names(Stunted.ent.2)<-c('Income','BMRMultiplier','Class')
RW.ent.all<-rbind(NARWnow.ent.2,NARW2000.ent.2,SRW.ent.2,Stunted.ent.2)
#Create a final category
RW.ent.all$Class2<-paste(RW.ent.all$BMRMultiplier,RW.ent.all$Class, sep="")
p16_r<-ggplot(RW.ent.all, aes(x=Income, y=Class2, color=Class2,
fill=Class2)) +
stat_halfeye(
point_color=NA, .width=0, height=0.6,
position=position_nudge(y = 0.3)
) +
geom_boxplot(
position=position_nudge(y = 0.2),
width=0.1, outlier.shape=NA, alpha=0.1
)+
xlim(0,40)+
labs(x=expression("Energy income as a multiplier of BMR (E)"))+
theme(axis.title.y= element_blank(), axis.text.y=element_blank(), axis.ticks.y=element_blank())+
scale_colour_manual(values = c("#F8766D","#F8766D","#F8766D","#F8766D","#7CAE00","#7CAE00","#7CAE00","#7CAE00"
,"#00BFC4","#00BFC4","#00BFC4","#00BFC4","#C77CFF","#C77CFF","#C77CFF","#C77CFF"))+
scale_fill_manual(values = c("#F8766D","#F8766D","#F8766D","#F8766D","#7CAE00","#7CAE00","#7CAE00","#7CAE00"
,"#00BFC4","#00BFC4","#00BFC4","#00BFC4","#C77CFF","#C77CFF","#C77CFF","#C77CFF"))+
scale_y_discrete(limits=rev,labels =
c("Stunted","SRW","NARWnow","NARW2000","Stunted","SRW","NARWnow","NARW2000","Stunted","SRW","NARWnow","NARW2000","
Stunted","SRW","NARWnow","NARW2000"))
p16_r
######################################################
#What do these values look like?
tapply(RW.ent.all$Income, RW.ent.all$Class2, summary)
RW.ent.all.classes <- RW.ent.all[RW.ent.all$Class=="NARWnow.ent"|RW.ent.all$Class=="NARW2000.ent"|RW.ent.all$Class=="SRW.ent",]
mean(RW.ent.all.classes$Income[RW.ent.all.classes$BMRMultiplier=="BMR*0.3"])
sd(RW.ent.all.classes$Income[RW.ent.all.classes$BMRMultiplier=="BMR*0.3"])
mean(RW.ent.all.classes$Income[RW.ent.all.classes$BMRMultiplier=="BMR*1"])
sd(RW.ent.all.classes$Income[RW.ent.all.classes$BMRMultiplier=="BMR*1"])
#reset the names of the food intake
names(Jin.ent2)<-c("kJintake","Class","MJintake")
######################################################
#Daily entanglement cost relative to food intake
Jin.Tngl.ent<-cbind(Jin.ent2,Tanglerange)
Jin.Tngl.ent$TangleToIntake<-Jin.Tngl.ent$Tanglerange/Jin.Tngl.ent$MJintake
#plot this
p17_r<-ggplot(Jin.Tngl.ent, aes(x=TangleToIntake, y=Class, color=Class,
fill=Class)) +
stat_halfeye(
point_color=NA, .width=0, height=0.6,
position=position_nudge(y = 0.3)
) +
geom_boxplot(
position=position_nudge(y = 0.2),
width=0.1, outlier.shape=NA, alpha=0.1
)+
xlim(0,4)+
labs(x=expression("Entanglement expenditure per energy income (E)"))+
scale_colour_manual(values = c("black", "black","black","black"))+
scale_fill_manual(values = c("black", "black","black","black"))+
theme(axis.title.y= element_blank(),axis.text.y=element_blank(), axis.ticks.y=element_blank())
p17_r
############################################
#These values as percentages for the three classes
#Estimates for only the three classes
Jin.Tngl.entClasses <- Jin.Tngl.ent[Jin.Tngl.ent$Class=="NARWnow"|Jin.Tngl.ent$Class=="NARW2000"|Jin.Tngl.ent$Class=="SRW",]
round(mean(Jin.Tngl.entClasses$TangleToIntake),2)*100
round(sd(Jin.Tngl.entClasses$TangleToIntake),3)*100
#95% CIs
Jin.Tngl.ent.sortedClasses<-Jin.Tngl.entClasses[order(Jin.Tngl.entClasses$TangleToIntake),]
#5%
round(Jin.Tngl.ent.sortedClasses[1500,5],3)*100
#95%
round(Jin.Tngl.ent.sortedClasses[28500,5],3)*100
#These values as percentages for all
round(mean(Jin.Tngl.ent$TangleToIntake),2)*100
round(sd(Jin.Tngl.ent$TangleToIntake),3)*100
#95% CIs
Jin.Tngl.ent.sorted<-Jin.Tngl.ent[order(Jin.Tngl.ent$TangleToIntake),]
#5%
round(Jin.Tngl.ent.sorted[2000,5],3)*100
#95%
round(Jin.Tngl.ent.sorted[38000,5],3)*100
#
round(min(Jin.Tngl.ent$TangleToIntake),3)*100
round(max(Jin.Tngl.ent$TangleToIntake),3)*100
aggregate(Jin.Tngl.ent$TangleToIntake, list(Jin.Tngl.ent$Class), FUN=mean)
aggregate(Jin.Tngl.ent$TangleToIntake, list(Jin.Tngl.ent$Class), FUN=sd)
#################################################
#Compare with standard feeding
round(mean(Jin.Tngl$TangleToIntake),2)
round(sd(Jin.Tngl$TangleToIntake),3)
aggregate(Jin.Tngl$TangleToIntake, list(Jin.Tngl$Class), FUN=mean)
aggregate(Jin.Tngl$TangleToIntake, list(Jin.Tngl$Class), FUN=sd)
############################################
############################################
#######Other comparative calculations#######
############################################
############################################
#sort data
comp <- as.data.frame(cbind(Jin2.sorted, Jin.ent2.sorted, Tanglerange))
names(comp)<-c("kJintakeu","Classu","MJintakeu","kJintakee","Classe","MJintakee","ET")
comp <- comp[-c(1,4,5)]
names(comp)<-c("Class","MJintakeu","MJintakee","ET")
##########Difference energy intake vs energy intake entangled
comp$intake.minus.intakee <- comp$MJintakeu-comp$MJintakee
mean(comp$intake.minus.intakee)
sd(comp$intake.minus.intakee)
comp.sorted.1 <- comp[order(comp$intake.minus.intakee),]
#5%
round(comp.sorted.1[2000,5],3)
#95%
round(comp.sorted.1[38000,5],3)
##########Energy intake entangled minus entanglement drag costs
comp$intakee.minus.ET <- comp$MJintakee-comp$ET
mean(comp$intakee.minus.ET)
sd(comp$intakee.minus.ET)
comp.sorted.2 <- comp[order(comp$intakee.minus.ET),]
#5%
round(comp.sorted.2[2000,6],3)
#95%
round(comp.sorted.2[38000,6],3)
##########Energy that can be spent when entangled in comparison to not entangled
comp$Edifference <- (comp$MJintakee-comp$ET)/comp$MJintakeu
mean(comp$Edifference)
sd(comp$Edifference)
comp.sorted.3 <- comp[order(comp$Edifference),]
#5%
round(comp.sorted.3[2000,7],3)
#95%
round(comp.sorted.3[38000,7],3)
################################################
################################################
########Final plot creation and saving##########
################################################
################################################
#Figure 1
Figure1 <- ggarrange(p1_r,p13_r,p3_r,p14_r,p4_r,p15_r,
ncol=2,
nrow=3)
Figure1
ggsave ("Output/Figure_1.pdf", Figure1,
width=16.9,
height=21,
units=c("cm"),
dpi = 1000)
ggsave ("Output/Figure_1.jpeg", Figure1,
width=16.9,
height=21,
units=c("cm"),
dpi = 1000)
#Figure 2
Figure2 <- ggarrange(p5,p12,
ncol=2, nrow=1, common.legend = TRUE)
Figure2
ggsave ("Output/Figure_2.pdf", Figure2,
width=16.9,
height=21,
units=c("cm"),
dpi = 1000)
ggsave ("Output/Figure_2.jpeg", Figure2,
width=16.9,
height=21,
units=c("cm"),
dpi = 1000)
#Figure 3
Figure3 <- ggarrange(p10,p16_r,p11,p17_r,
ncol=2,
nrow=2)
Figure3
ggsave ("Output/Figure_3.pdf", Figure3,
width=16.9,
height=28,
units=c("cm"),
dpi = 1000)
ggsave ("Output/Figure_3.jpeg", Figure3,
width=16.9,
height=28,
units=c("cm"),
dpi = 1000)
##################################################
##################References######################
##################################################
#Baumgartner MF, Mate BR (2003) Summertime foraging ecology of North Atlantic right whales. Marine Ecology Progress Series 264:123-135
#Christiansen F, Dawson SM, Durban JW, Fearnbach H, Miller CA, Bejder L, Uhart M, Sironi M, Corkeron P, Rayment W, Leunissen E, Haria E,
Ward R, Warick HA, Kerr I, Lynn MS, Pettis HM, Moore MJ (2020) Population comparison of right whale body condition reveals poor state of the
North Atlantic right whale. Marine Ecology Progress Series 640:1-16
#Kenney RD, Hyman MAM, Owen RE, Scott GP, Winn HE (1986) Estimation of Prey Densities Required by Western North Atlantic Right
Whales. Marine Mammal Science 2:1-13
#Kleiber M (1975) The fire of life: an introduction to animal energetics. R. E. Krieger Pub. Co., Huntington, N.Y.
#Maps F, Plourde S, Zakardjian B (2010) Control of dormancy by lipid metabolism in Calanus finmarchicus: a population model test. Marine
Ecology Progress Series 403:165-180
#Maps F, Zakardjian BA, Plourde S, Saucier FJ (2011) Modeling the interactions between the seasonal and diel migration behaviors of Calanus
finmarchicus and the circulation in the Gulf of St. Lawrence (Canada). Journal of Marine Systems 88:183-202
#McKinstry CAE, Westgate AJ, Koopman HN (2013) Annual variation in the nutritional value of Stage V Calanus finmarchicus: implications for
right whales and other copepod predators. Endangered Species Research 20:195-204
#Michaud J, Taggart CT (2007) Lipid and gross energy content of North Atlantic right whale food, Calanus finmarchicus, in the Bay of Fundy.
Endangered Species Research 3:77-94
#Miller CA, Best PB, Perryman WL, Baumgartner MF, Moore MJ (2012) Body shape changes associated with reproductive status, nutritive
condition and growth in right whales Eubalaena glacialis and E. australis. Marine Ecology Progress Series 459:135-156
#Stewart JD, Durban JW, Knowlton AR, Lynn MS, Fearnbach H, Barbaro J, Perryman WL, Miller CA, Moore MJ (2021) Decreasing body lengths
in North Atlantic right whales. Current Biology
#van der Hoop JM, Corkeron PJ, Moore MJ (2017) Entanglement is a costly life-history stage in large whales. Ecology and Evolution 7:92-106
#Wiebe PH, Boyd SH, Cox JL (1988) Functional regression equations for zooplankton displacement volume, wet weight, dry weight, and carbon: a
correction. Fishery Bulletin 86:833-835