Ecological-lifestyles-and-the-scaling-of-shark-gll-surface-area/Bigman et al 2018 J Morph Code.R at master · jennybigman/Ecological-lifestyles-and-the-scaling-of-shark-gll-surface-area · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
# Code for analyses in Bigman et al. (2018) Journal of Morphology

# Please cite the following if code or data is used:
# Bigman JS, Pardo SA, Prinzing TS, Dando M, Wegner NC, Dulvy NK.
# Ecological lifestyles and the scaling of shark gill surface area.
# Journal of Morphology. 2018;1–9. https://doi.org/10.1002/jmor.20879

# packages
# devtools::install_github("dgrtwo/broom")
library(tidyverse)
library(tidyr)
library(cowplot)
library(ggplot2)
library(MASS)
library(broom)
library(lme4)
library(nlme)
library(forcats)

#### Data wrangling ####

## Mustelus gill surface area data
MusCaliGSA <- read.csv("MustelusCalifornicusGSA_Data.csv", header = TRUE,
                       strip.white = TRUE) %>% tbl_df() # smoothhound raw data
## Gill surface area data for other 11 sharks
AllSharkData <- read.csv("Bigman_raw_data_combined.csv", header = TRUE,
                         strip.white = TRUE) %>% tbl_df

## Mustelus data transformations

#Gill surface area
MusCaliGSA$MassG <- (MusCaliGSA$MassKG*1000)
MusCaliGSA$Log10MassG <- log10(MusCaliGSA$MassG)
MusCaliGSA$Log10GSAcm2 <- log10(MusCaliGSA$GSAcm2)

# Total filament length
MusCaliGSA$Log10FilamentLength <- log10(MusCaliGSA$TotalFilamentLength2Sides)

# Average lamellar frequency
MusCaliGSA$Log10LamFreq <- log10(MusCaliGSA$LamellarFrequency)

# Mean bilateral lamellar surface area
MusCaliGSA$Log10LamSA <- log10(MusCaliGSA$BilateralLamellarSA)


## Remaining shark data transformations
AllSharkData$Log10MassG <- log10(AllSharkData$MassG)
AllSharkData$Log10GSAcm2 <- log10(AllSharkData$GSAcm2)

AllSharkData$Log10CenterMassG <- scale(AllSharkData$Log10MassG,
                                       center = log10(5000), scale = FALSE)


#### Analyses ####

## Mustelus gill surface area

# Total gill surface area regression and prediction
GSA.mass.log <- lm(Log10GSAcm2 ~ Log10MassG, data = MusCaliGSA)
summary(GSA.mass.log)

predict <- predict(GSA.mass.log, interval="prediction")
MusCaliGSA <- cbind(MusCaliGSA, predict)
MusCaliGSA$Fit.NoLog <- 10^(MusCaliGSA$fit)
MusCaliGSA$Upper.NoLog <- 10^(MusCaliGSA$upr)
MusCaliGSA$Lower.NoLog <- 10^(MusCaliGSA$lwr)


# Total filament length regression and prediction
FilamentLength.mass <- lm(Log10FilamentLength ~ Log10MassG, data = MusCaliGSA)
summary(FilamentLength.mass)

predict.FilamentLengthLog10 <- predict(FilamentLength.mass, interval="prediction")
MusCaliGSA.FilamentLength <- data.frame(MusCaliGSA$MassG, MusCaliGSA$TotalFilamentLength2Sides,
                                        MusCaliGSA$Log10FilamentLength, MusCaliGSA$Log10MassG,
                                        predict.FilamentLengthLog10)
MusCaliGSA.FilamentLength$FitFilLength.NoLog <- 10^(MusCaliGSA.FilamentLength$fit)
MusCaliGSA.FilamentLength$UpperFilLength.NoLog <- 10^(MusCaliGSA.FilamentLength$upr)
MusCaliGSA.FilamentLength$LowerFilLength.NoLog <- 10^(MusCaliGSA.FilamentLength$lwr)


# Average lamellar frequency regression and prediction
LamFreq.mass <- lm(Log10LamFreq ~ Log10MassG, data = MusCaliGSA)
summary(LamFreq.mass)

predict.Log10LamFreq <- predict(LamFreq.mass, interval="prediction")
MusCaliGSA.LamFreq <-data.frame(MusCaliGSA$MassG, MusCaliGSA$LamellarFrequency,
                                MusCaliGSA$Log10MassG, MusCaliGSA$Log10LamFreq,
                                predict.Log10LamFreq)
MusCaliGSA.LamFreq$FitFilLength.NoLog <- 10^(MusCaliGSA.LamFreq$fit)
MusCaliGSA.LamFreq$UpperFilLength.NoLog <- 10^(MusCaliGSA.LamFreq$upr)
MusCaliGSA.LamFreq$LowerFilLength.NoLog <- 10^(MusCaliGSA.LamFreq$lwr)


# Mean bilateral lamellar surface area regression and prediction
LamSA.mass <- lm(Log10LamSA ~ Log10MassG, data = MusCaliGSA)
summary(LamSA.mass)

predict.Log10LamSA <- predict(LamSA.mass, interval="prediction")
MusCaliGSA.LamSA <- data.frame(MusCaliGSA$MassG, MusCaliGSA$BilateralLamellarSA,
                               MusCaliGSA$Log10LamSA, MusCaliGSA$Log10MassG,
                               predict.Log10LamSA)
MusCaliGSA.LamSA$FitFilLength.NoLog <- 10^(MusCaliGSA.LamSA$fit)
MusCaliGSA.LamSA$UpperFilLength.NoLog <- 10^(MusCaliGSA.LamSA$upr)
MusCaliGSA.LamSA$LowerFilLength.NoLog <- 10^(MusCaliGSA.LamSA$lwr)


#### Scaling of gill surface area for all 12 sharks ####

## Regressions to estimate coefficients
fitted_models.Log10CenterMass <- AllSharkData %>%
  group_by(Species) %>%
  do(fits.Log10CenterMass = lm(Log10GSAcm2 ~ Log10CenterMassG, data = .))

model.output.Log10CenterMass <- fitted_models.Log10CenterMass %>%
  tidy(fits.Log10CenterMass)
fits.Log10CenterMass <- fitted_models.Log10CenterMass %>%
  augment(fits.Log10CenterMass)

## Do intercepts and slopes differ across species?
mod1 <- lm(Log10GSAcm2 ~ Log10CenterMassG * Species - 1, data = AllSharkData) #
coef(mod1)
summary(mod1)

## Do ecological lifestyle traits differ across species?

# 1. Caudal fin aspect ratio
ActivityLevelDifferences  <- lme(Log10GSAcm2 ~ Log10CenterMassG * AspectRatio,
                                 random = ~ Log10CenterMassG | Species,
                                 data = AllSharkData, control = lmeControl(opt = "optim"))
coef(ActivityLevelDifferences)
summary(ActivityLevelDifferences)

# 2. Habitat type
HabitatDifferences  <- lme(Log10GSAcm2 ~ Log10CenterMassG * HabitatType,
                           random = ~ Log10CenterMassG | Species,
                           data = AllSharkData, control = lmeControl(opt = "optim"))
coef(HabitatDifferences)
summary(HabitatDifferences)

# 3. Maximum size
MaxSizeDifferences  <- lme(Log10GSAcm2 ~ Log10CenterMassG * LogMaxMassG,
                           random = ~ Log10CenterMassG | Species,
                           data = AllSharkData, control = lmeControl(opt = "optim"))
coef(MaxSizeDifferences)
summary(MaxSizeDifferences)