L. buccanella body size models
Naomi Bakken
2024-08-27
This R markdown file includes all the R code used to analyze Blackfin Snapper (Lutjanus buccanella) fork length-weight, total length-weight, and maturity relationships.
Load libraries
library(tidyverse)
library(brms) # for Bayesian regression modeling using Stan
library(rstudioapi)
library(nls.multstart) # for nonlinear least squares regression modeling with multiple start values
library(stringr) # for wrapping text on graphics
library(vtable)
library(viridis)Define color palettes
safe_colorblind_palette <- c("#88CCEE", "#CC6677", "#DDCC77", "#117733", "#332288", "#AA4499", "#44AA99", "#999933", "#882255", "#661100", "#6699CC", "#888888")
turbo_colors_onboard <- viridis_pal(option = "turbo", begin = 0.05, end = 0.95)(9)
turbo_colors_onboard## [1] "#3E378FFF" "#4683F9FF" "#1FC9DDFF" "#3BF58FFF" "#A2FC3CFF" "#E7D739FF"
## [7] "#FE932AFF" "#E3440AFF" "#A31301FF"Load data file
bsdf_2506_LLrelationships <- read.csv("bs2506df.csv")
bsdf_2506_LLrelationships$catch_month <- factor(bsdf_2506_LLrelationships$catch_month)
bsdf_2506_LLrelationships$year <- factor(bsdf_2506_LLrelationships$year)
bsdf_2506_LLrelationships$species <- factor(bsdf_2506_LLrelationships$species)
bsdf_2506_LLrelationships$common_name <- factor(bsdf_2506_LLrelationships$common_name)
bsdf_2506_LLrelationships$frozen <- factor(bsdf_2506_LLrelationships$frozen)
bsdf_2506_LLrelationships$catch_method <- factor(bsdf_2506_LLrelationships$catch_method)
bsdf_2506_LLrelationships$bait <- factor(bsdf_2506_LLrelationships$bait)
bsdf_2506_LLrelationships$area <- factor(bsdf_2506_LLrelationships$area)
bsdf_2506_LLrelationships$gutted <- factor(bsdf_2506_LLrelationships$gutted)
bsdf_2506_LLrelationships$sex <- factor(bsdf_2506_LLrelationships$sex)
bsdf_2506_LLrelationships$maturity <- factor(bsdf_2506_LLrelationships$maturity)
summary(bsdf_2506_LLrelationships)## X catch_date catch_month landing_time
## Min. : 1 Length:85 : 3 Length:85
## 1st Qu.:22 Class :character Apr:11 Class :character
## Median :43 Mode :character Dec: 1 Mode :character
## Mean :43 Jul: 9
## 3rd Qu.:64 Jun:32
## Max. :85 May:19
## Sep:10
## dissection_date dissection_time year month
## Length:85 Length:85 2023:52 Length:85
## Class :character Class :character 2024:33 Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
## observer species common_name
## Length:85 Lutjanus buccanella:85 Blackfin Snapper:85
## Class :character
## Mode :character
##
##
##
##
## fish_code frozen catch_method bait
## Length:85 N:33 :51 :51
## Class :character Y:52 Lobster trap:11 Cowskin :11
## Mode :character Redfish trap:23 Japanese bait :22
## Japanese bait; Squid: 1
##
##
##
## area SL_CM FL_CM TL_CM HH_CM
## :41 Min. :16.40 Min. :18.70 Min. :20.50 Min. :3.600
## A3/A4 :17 1st Qu.:24.60 1st Qu.:23.60 1st Qu.:26.00 1st Qu.:6.200
## C2 :10 Median :27.95 Median :26.30 Median :28.40 Median :7.950
## D3 :10 Mean :28.33 Mean :27.66 Mean :29.89 Mean :7.359
## D5 : 3 3rd Qu.:34.08 3rd Qu.:30.00 3rd Qu.:32.48 3rd Qu.:8.850
## D4 : 2 Max. :38.30 Max. :42.90 Max. :45.90 Max. :9.900
## (Other): 2 NA's :51 NA's :3 NA's :63
## ED_CM W_G finclip otoliths gutted
## Min. :1.200 Min. : 108.5 Length:85 Length:85 N:84
## 1st Qu.:2.200 1st Qu.: 234.7 Class :character Class :character Y: 1
## Median :2.400 Median : 322.9 Mode :character Mode :character
## Mean :2.245 Mean : 428.9
## 3rd Qu.:2.500 3rd Qu.: 471.1
## Max. :2.900 Max. :1493.9
## NA's :63 NA's :3
## gonads_present sex maturity G1L_CM G2L_CM
## Length:85 :32 :52 Min. :1.100 Min. :1.200
## Class :character F:19 A:15 1st Qu.:2.525 1st Qu.:2.575
## Mode :character M:34 B:10 Median :3.200 Median :3.200
## C: 8 Mean :3.698 Mean :3.720
## 3rd Qu.:4.200 3rd Qu.:4.525
## Max. :9.300 Max. :8.800
## NA's :19 NA's :25
## GW_G stomach._sample
## Length:85 Length:85
## Class :character Class :character
## Mode :character Mode :character
##
##
##
## View(bsdf_2506_LLrelationships)
sumtable(bsdf_2506_LLrelationships, digits = 4)| Variable | N | Mean | Std. Dev. | Min | Pctl. 25 | Pctl. 75 | Max |
|---|---|---|---|---|---|---|---|
| X | 85 | 43 | 24.68 | 1 | 22 | 64 | 85 |
| catch_month | 85 | ||||||
| … | 3 | 3.53% | |||||
| … Apr | 11 | 12.94% | |||||
| … Dec | 1 | 1.18% | |||||
| … Jul | 9 | 10.59% | |||||
| … Jun | 32 | 37.65% | |||||
| … May | 19 | 22.35% | |||||
| … Sep | 10 | 11.76% | |||||
| year | 85 | ||||||
| … 2023 | 52 | 61.18% | |||||
| … 2024 | 33 | 38.82% | |||||
| observer | 85 | ||||||
| … HK | 37 | 43.53% | |||||
| … LJ | 8 | 9.41% | |||||
| … NB | 9 | 10.59% | |||||
| … QH | 7 | 8.24% | |||||
| … SM | 15 | 17.65% | |||||
| … TO | 9 | 10.59% | |||||
| species | 85 | ||||||
| … Lutjanus buccanella | 85 | 100% | |||||
| common_name | 85 | ||||||
| … Blackfin Snapper | 85 | 100% | |||||
| frozen | 85 | ||||||
| … N | 33 | 38.82% | |||||
| … Y | 52 | 61.18% | |||||
| catch_method | 85 | ||||||
| … | 51 | 60% | |||||
| … Lobster trap | 11 | 12.94% | |||||
| … Redfish trap | 23 | 27.06% | |||||
| bait | 85 | ||||||
| … | 51 | 60% | |||||
| … Cowskin | 11 | 12.94% | |||||
| … Japanese bait | 22 | 25.88% | |||||
| … Japanese bait; Squid | 1 | 1.18% | |||||
| area | 85 | ||||||
| … | 41 | 48.24% | |||||
| … A3/A4 | 17 | 20% | |||||
| … A4/A5 | 1 | 1.18% | |||||
| … B4 | 1 | 1.18% | |||||
| … C2 | 10 | 11.76% | |||||
| … D3 | 10 | 11.76% | |||||
| … D4 | 2 | 2.35% | |||||
| … D5 | 3 | 3.53% | |||||
| SL_CM | 34 | 28.33 | 5.769 | 16.4 | 24.6 | 34.08 | 38.3 |
| FL_CM | 85 | 27.66 | 5.853 | 18.7 | 23.6 | 30 | 42.9 |
| TL_CM | 82 | 29.89 | 6.296 | 20.5 | 26 | 32.48 | 45.9 |
| HH_CM | 22 | 7.359 | 1.929 | 3.6 | 6.2 | 8.85 | 9.9 |
| ED_CM | 22 | 2.245 | 0.424 | 1.2 | 2.2 | 2.5 | 2.9 |
| W_G | 82 | 428.9 | 321.1 | 108.5 | 234.7 | 471.1 | 1494 |
| finclip | 85 | ||||||
| … | 49 | 57.65% | |||||
| … N | 10 | 11.76% | |||||
| … Y | 26 | 30.59% | |||||
| otoliths | 85 | ||||||
| … N | 3 | 3.53% | |||||
| … y | 1 | 1.18% | |||||
| … Y | 81 | 95.29% | |||||
| gutted | 85 | ||||||
| … N | 84 | 98.82% | |||||
| … Y | 1 | 1.18% | |||||
| gonads_present | 85 | ||||||
| … N | 18 | 21.18% | |||||
| … Y | 67 | 78.82% | |||||
| sex | 85 | ||||||
| … | 32 | 37.65% | |||||
| … F | 19 | 22.35% | |||||
| … M | 34 | 40% | |||||
| maturity | 85 | ||||||
| … | 52 | 61.18% | |||||
| … A | 15 | 17.65% | |||||
| … B | 10 | 11.76% | |||||
| … C | 8 | 9.41% | |||||
| G1L_CM | 66 | 3.698 | 1.761 | 1.1 | 2.525 | 4.2 | 9.3 |
| G2L_CM | 60 | 3.72 | 1.713 | 1.2 | 2.575 | 4.525 | 8.8 |
| stomach._sample | 85 | ||||||
| … N | 70 | 82.35% | |||||
| … Y | 15 | 17.65% |
Fork length (FL_CM) - Total length (TL_CM) relationship
bs2506_FLTL_plotbase <- ggplot(bsdf_2506_LLrelationships, aes(x = FL_CM, y = TL_CM)) +
geom_point(aes(fill = sex),
pch = 21,
size = 2) +
theme_classic() +
labs(x = "Fork Length in cm",
y = "Total Length in cm",
title = "Lutjanus buccanella FL-TL relationship (n = 82)",
fill = "Sex") +
scale_fill_viridis(option = "C",
discrete = TRUE,
begin = 0.2,
end = 0.8,
labels = c("N/A", "Female", "Male"))
bs2506_FLTL_plotbase
Linear model
bs2506_FLTL_lm <- lm(FL_CM ~ TL_CM, data = bsdf_2506_LLrelationships)
summary(bs2506_FLTL_lm)##
## Call:
## lm(formula = FL_CM ~ TL_CM, data = bsdf_2506_LLrelationships)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.5541 -0.1767 -0.0475 0.2763 1.8519
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.43743 0.37181 -1.176 0.243
## TL_CM 0.94021 0.01218 77.211 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.69 on 80 degrees of freedom
## (3 observations deleted due to missingness)
## Multiple R-squared: 0.9868, Adjusted R-squared: 0.9866
## F-statistic: 5962 on 1 and 80 DF, p-value: < 2.2e-16length(bs2506_FLTL_lm$fitted.values)## [1] 82Scatterplot the residuals against the fitted values to check homogeneity of variance
plot(resid(bs2506_FLTL_lm) ~ fitted(bs2506_FLTL_lm))
Fine if you disregard the one outlier
Plot a histogram to check normality of residuals
qqnorm(bs2506_FLTL_lm$residuals)
qqline(bs2506_FLTL_lm$residuals)
Fork length (FL_CM) - Standard length (SL_CM) relationship
bs2506_FLSL_plot <- ggplot(bsdf_2506_LLrelationships, aes(x = FL_CM, y = SL_CM)) +
geom_point(aes(fill = sex),
pch = 21,
size = 2) +
theme_classic() +
labs(x = "Fork Length in cm",
y = "Standard Length in cm",
title = "Lutjanus buccanella FL-SL relationship (n = 34)",
fill = "Sex") +
scale_fill_viridis(option = "C",
discrete = TRUE,
begin = 0.2,
end = 0.8,
labels = c("N/A", "Female", "Male"))
bs2506_FLSL_plot
Linear model
bs2506_FLSL_lm <- lm(FL_CM ~ SL_CM, data = bsdf_2506_LLrelationships)
summary(bs2506_FLSL_lm)##
## Call:
## lm(formula = FL_CM ~ SL_CM, data = bsdf_2506_LLrelationships)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.04176 -0.39552 0.07946 0.28097 0.92115
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.14872 0.44344 0.335 0.74
## SL_CM 1.12878 0.01535 73.550 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5086 on 32 degrees of freedom
## (51 observations deleted due to missingness)
## Multiple R-squared: 0.9941, Adjusted R-squared: 0.9939
## F-statistic: 5410 on 1 and 32 DF, p-value: < 2.2e-16length(bs2506_FLSL_lm$fitted.values)## [1] 34Scatterplot the residuals against the fitted values to check homogeneity of variance
plot(resid(bs2506_FLSL_lm) ~ fitted(bs2506_FLSL_lm))
Fine
Plot a histogram to check normality of residuals
qqnorm(bs2506_FLSL_lm$residuals)
qqline(bs2506_FLSL_lm$residuals)
Fork length (TL_CM) - Standard length (SL_CM) relationship
bs2506_TLSL_plot <- ggplot(bsdf_2506_LLrelationships, aes(x = TL_CM, y = SL_CM)) +
geom_point(aes(fill = sex),
pch = 21,
size = 2) +
theme_classic() +
labs(x = "Total Length in cm",
y = "Standard Length in cm",
title = "Lutjanus buccanella TL-SL relationship (n = 34)",
fill = "Sex") +
scale_fill_viridis(option = "C",
discrete = TRUE,
begin = 0.2,
end = 0.8,
labels = c("N/A", "Female", "Male"))
bs2506_FLSL_plot
Linear model
bs2506_TLSL_lm <- lm(TL_CM ~ SL_CM, data = bsdf_2506_LLrelationships)
summary(bs2506_TLSL_lm)##
## Call:
## lm(formula = TL_CM ~ SL_CM, data = bsdf_2506_LLrelationships)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.20835 -0.50237 0.02702 0.41061 1.13653
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.44740 0.53540 0.836 0.41
## SL_CM 1.20244 0.01853 64.892 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6141 on 32 degrees of freedom
## (51 observations deleted due to missingness)
## Multiple R-squared: 0.9925, Adjusted R-squared: 0.9922
## F-statistic: 4211 on 1 and 32 DF, p-value: < 2.2e-16length(bs2506_TLSL_lm$fitted.values)## [1] 34Scatterplot the residuals against the fitted values to check homogeneity of variance
plot(resid(bs2506_TLSL_lm) ~ fitted(bs2506_TLSL_lm))
Fine
Plot a histogram to check normality of residuals
qqnorm(bs2506_TLSL_lm$residuals)
qqline(bs2506_TLSL_lm$residuals) 
Fork length (FL_CM) - Weight (W_G) relationship
Load data file
bsdf_2506_FL_W <- read.csv("bs2506df.csv")
bsdf_2506_FL_W$catch_month <- factor(bsdf_2506_FL_W$catch_month)
bsdf_2506_FL_W$year <- factor(bsdf_2506_FL_W$year)
bsdf_2506_FL_W$species <- factor(bsdf_2506_FL_W$species)
bsdf_2506_FL_W$common_name <- factor(bsdf_2506_FL_W$common_name)
bsdf_2506_FL_W$frozen <- factor(bsdf_2506_FL_W$frozen)
bsdf_2506_FL_W$catch_method <- factor(bsdf_2506_FL_W$catch_method)
bsdf_2506_FL_W$bait <- factor(bsdf_2506_FL_W$bait)
bsdf_2506_FL_W$area <- factor(bsdf_2506_FL_W$area)
bsdf_2506_FL_W$gutted <- factor(bsdf_2506_FL_W$gutted)
bsdf_2506_FL_W$sex <- factor(bsdf_2506_FL_W$sex)
bsdf_2506_FL_W$maturity <- factor(bsdf_2506_FL_W$maturity)
# Remove gutted samples
bsdf_2506_FL_W <- subset(filter(bsdf_2506_FL_W, gutted == "N"))
View(bsdf_2506_FL_W)
summary(bsdf_2506_FL_W)## X catch_date catch_month landing_time
## Min. : 2.00 Length:84 : 3 Length:84
## 1st Qu.:22.75 Class :character Apr:10 Class :character
## Median :43.50 Mode :character Dec: 1 Mode :character
## Mean :43.50 Jul: 9
## 3rd Qu.:64.25 Jun:32
## Max. :85.00 May:19
## Sep:10
## dissection_date dissection_time year month
## Length:84 Length:84 2023:52 Length:84
## Class :character Class :character 2024:32 Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
## observer species common_name
## Length:84 Lutjanus buccanella:84 Blackfin Snapper:84
## Class :character
## Mode :character
##
##
##
##
## fish_code frozen catch_method bait
## Length:84 N:32 :51 :51
## Class :character Y:52 Lobster trap:10 Cowskin :10
## Mode :character Redfish trap:23 Japanese bait :22
## Japanese bait; Squid: 1
##
##
##
## area SL_CM FL_CM TL_CM HH_CM
## :41 Min. :16.40 Min. :18.70 Min. :20.50 Min. :3.600
## A3/A4 :17 1st Qu.:24.60 1st Qu.:23.60 1st Qu.:26.00 1st Qu.:6.200
## C2 :10 Median :28.00 Median :26.25 Median :28.30 Median :7.950
## D3 :10 Mean :28.44 Mean :27.67 Mean :29.89 Mean :7.359
## D5 : 3 3rd Qu.:34.20 3rd Qu.:30.05 3rd Qu.:32.50 3rd Qu.:8.850
## D4 : 2 Max. :38.30 Max. :42.90 Max. :45.90 Max. :9.900
## (Other): 1 NA's :51 NA's :3 NA's :62
## ED_CM W_G finclip otoliths gutted
## Min. :1.200 Min. : 108.5 Length:84 Length:84 N:84
## 1st Qu.:2.200 1st Qu.: 234.3 Class :character Class :character Y: 0
## Median :2.400 Median : 323.4 Mode :character Mode :character
## Mean :2.245 Mean : 430.3
## 3rd Qu.:2.500 3rd Qu.: 474.0
## Max. :2.900 Max. :1493.9
## NA's :62 NA's :3
## gonads_present sex maturity G1L_CM G2L_CM
## Length:84 :32 :52 Min. :1.100 Min. :1.200
## Class :character F:18 A:15 1st Qu.:2.500 1st Qu.:2.550
## Mode :character M:34 B:10 Median :3.200 Median :3.200
## C: 7 Mean :3.709 Mean :3.705
## 3rd Qu.:4.200 3rd Qu.:4.500
## Max. :9.300 Max. :8.800
## NA's :19 NA's :25
## GW_G stomach._sample
## Length:84 Length:84
## Class :character Class :character
## Mode :character Mode :character
##
##
##
## # remove NA rows
bsdf_2506_FL_W <- bsdf_2506_FL_W%>%
drop_na(FL_CM)
View(bsdf_2506_FL_W)
summary(bsdf_2506_FL_W)## X catch_date catch_month landing_time
## Min. : 2.00 Length:84 : 3 Length:84
## 1st Qu.:22.75 Class :character Apr:10 Class :character
## Median :43.50 Mode :character Dec: 1 Mode :character
## Mean :43.50 Jul: 9
## 3rd Qu.:64.25 Jun:32
## Max. :85.00 May:19
## Sep:10
## dissection_date dissection_time year month
## Length:84 Length:84 2023:52 Length:84
## Class :character Class :character 2024:32 Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
## observer species common_name
## Length:84 Lutjanus buccanella:84 Blackfin Snapper:84
## Class :character
## Mode :character
##
##
##
##
## fish_code frozen catch_method bait
## Length:84 N:32 :51 :51
## Class :character Y:52 Lobster trap:10 Cowskin :10
## Mode :character Redfish trap:23 Japanese bait :22
## Japanese bait; Squid: 1
##
##
##
## area SL_CM FL_CM TL_CM HH_CM
## :41 Min. :16.40 Min. :18.70 Min. :20.50 Min. :3.600
## A3/A4 :17 1st Qu.:24.60 1st Qu.:23.60 1st Qu.:26.00 1st Qu.:6.200
## C2 :10 Median :28.00 Median :26.25 Median :28.30 Median :7.950
## D3 :10 Mean :28.44 Mean :27.67 Mean :29.89 Mean :7.359
## D5 : 3 3rd Qu.:34.20 3rd Qu.:30.05 3rd Qu.:32.50 3rd Qu.:8.850
## D4 : 2 Max. :38.30 Max. :42.90 Max. :45.90 Max. :9.900
## (Other): 1 NA's :51 NA's :3 NA's :62
## ED_CM W_G finclip otoliths gutted
## Min. :1.200 Min. : 108.5 Length:84 Length:84 N:84
## 1st Qu.:2.200 1st Qu.: 234.3 Class :character Class :character Y: 0
## Median :2.400 Median : 323.4 Mode :character Mode :character
## Mean :2.245 Mean : 430.3
## 3rd Qu.:2.500 3rd Qu.: 474.0
## Max. :2.900 Max. :1493.9
## NA's :62 NA's :3
## gonads_present sex maturity G1L_CM G2L_CM
## Length:84 :32 :52 Min. :1.100 Min. :1.200
## Class :character F:18 A:15 1st Qu.:2.500 1st Qu.:2.550
## Mode :character M:34 B:10 Median :3.200 Median :3.200
## C: 7 Mean :3.709 Mean :3.705
## 3rd Qu.:4.200 3rd Qu.:4.500
## Max. :9.300 Max. :8.800
## NA's :19 NA's :25
## GW_G stomach._sample
## Length:84 Length:84
## Class :character Class :character
## Mode :character Mode :character
##
##
##
## Remove gutted samples
bsdf_2506_FL_W <- subset(filter(bsdf_2506_FL_W, gutted == "N"))
View(bsdf_2506_FL_W)
summary(bsdf_2506_FL_W)## X catch_date catch_month landing_time
## Min. : 2.00 Length:84 : 3 Length:84
## 1st Qu.:22.75 Class :character Apr:10 Class :character
## Median :43.50 Mode :character Dec: 1 Mode :character
## Mean :43.50 Jul: 9
## 3rd Qu.:64.25 Jun:32
## Max. :85.00 May:19
## Sep:10
## dissection_date dissection_time year month
## Length:84 Length:84 2023:52 Length:84
## Class :character Class :character 2024:32 Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
## observer species common_name
## Length:84 Lutjanus buccanella:84 Blackfin Snapper:84
## Class :character
## Mode :character
##
##
##
##
## fish_code frozen catch_method bait
## Length:84 N:32 :51 :51
## Class :character Y:52 Lobster trap:10 Cowskin :10
## Mode :character Redfish trap:23 Japanese bait :22
## Japanese bait; Squid: 1
##
##
##
## area SL_CM FL_CM TL_CM HH_CM
## :41 Min. :16.40 Min. :18.70 Min. :20.50 Min. :3.600
## A3/A4 :17 1st Qu.:24.60 1st Qu.:23.60 1st Qu.:26.00 1st Qu.:6.200
## C2 :10 Median :28.00 Median :26.25 Median :28.30 Median :7.950
## D3 :10 Mean :28.44 Mean :27.67 Mean :29.89 Mean :7.359
## D5 : 3 3rd Qu.:34.20 3rd Qu.:30.05 3rd Qu.:32.50 3rd Qu.:8.850
## D4 : 2 Max. :38.30 Max. :42.90 Max. :45.90 Max. :9.900
## (Other): 1 NA's :51 NA's :3 NA's :62
## ED_CM W_G finclip otoliths gutted
## Min. :1.200 Min. : 108.5 Length:84 Length:84 N:84
## 1st Qu.:2.200 1st Qu.: 234.3 Class :character Class :character Y: 0
## Median :2.400 Median : 323.4 Mode :character Mode :character
## Mean :2.245 Mean : 430.3
## 3rd Qu.:2.500 3rd Qu.: 474.0
## Max. :2.900 Max. :1493.9
## NA's :62 NA's :3
## gonads_present sex maturity G1L_CM G2L_CM
## Length:84 :32 :52 Min. :1.100 Min. :1.200
## Class :character F:18 A:15 1st Qu.:2.500 1st Qu.:2.550
## Mode :character M:34 B:10 Median :3.200 Median :3.200
## C: 7 Mean :3.709 Mean :3.705
## 3rd Qu.:4.200 3rd Qu.:4.500
## Max. :9.300 Max. :8.800
## NA's :19 NA's :25
## GW_G stomach._sample
## Length:84 Length:84
## Class :character Class :character
## Mode :character Mode :character
##
##
##
## First look at the data with a scatterplot
fill = sex
ggplot(bsdf_2506_FL_W, aes(x = FL_CM, y = W_G)) +
geom_point(aes(fill = sex), color = "#1B1B1B", pch = 21, size = 2) +
theme_classic() +
labs(x = "Fork Length in cm",
y = "Weight in grams",
title = "Lutjanus buccanella, n = 81") +
scale_fill_viridis(option = "D", discrete = TRUE)
Curves from published literature on L. buccanella (BS) TL-W relationships
W_kg values calculated from Burton et al. (2016) parameters then converted to cm and g
Use size range from their paper (TL = 180mm-609mm) W = 9.54e-9FL^(3.11)
#predicted values
bs_burtonFL_x <- c(180:609)
bs_burtonFL_df <- data.frame(bs_burtonFL_x)
bs_burtonFL_f <- function(bs_burtonFL_x) {
# Replace this with your actual formula
return(9.54e-9*bs_burtonFL_x^3.11)
}
bs_burtonFL_df$bs_burtonFL_y <- sapply(bs_burtonFL_df$bs_burtonFL_x, bs_burtonFL_f)
#change kg to g
bs_burtonFL_df$bs_burtonFL_y_g <- bs_burtonFL_df$bs_burtonFL_y*1000
#change mm to cm
bs_burtonFL_df$bs_burtonFL_x_cm <- bs_burtonFL_df$bs_burtonFL_x/10
View(bs_burtonFL_df)Plot models from literature
bs_litFL_plottemp <- ggplot(bsdf_2506_FL_W, aes(x = FL_CM, y = W_G)) +
geom_point(aes(TL_CM, W_G),
color = "#242124",
alpha = 0.8) +
geom_smooth(aes(x = bs_burtonFL_x_cm, y = bs_burtonFL_y_g),
data = bs_burtonFL_df,
method = "loess",
se = FALSE,
linewidth = 0.75,
linetype = 2,
color = safe_colorblind_palette[1]) +
labs(x = "Fork Length in cm",
y = "Weight in grams",
title = "Blackfin Snapper Fork Length: Saba vs. Other Regions (n = 81)") +
theme_classic()
# add label with equation 7.79e-9*bs_burton_x^3.09
bs_litFL_plot <- bs_litFL_plottemp +
geom_text(x = 18, y = 2800,
label = expression(W == 17.79 %*% 10^{-09} * TL^{3.09} ~ "(SE & Caribbean U.S.)"),
colour = safe_colorblind_palette[1],
hjust = 0)
bs_litFL_plot
Linear model
bs2506_FL_lm1 <- lm(FL_CM ~ W_G, data = bsdf_2506_FL_W)
bs2506_FL_lm1##
## Call:
## lm(formula = FL_CM ~ W_G, data = bsdf_2506_FL_W)
##
## Coefficients:
## (Intercept) W_G
## 19.92836 0.01798summary(bs2506_FL_lm1)##
## Call:
## lm(formula = FL_CM ~ W_G, data = bsdf_2506_FL_W)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.3942 -0.8990 0.1493 0.9910 2.9111
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.993e+01 2.805e-01 71.03 <2e-16 ***
## W_G 1.798e-02 5.227e-04 34.41 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.509 on 79 degrees of freedom
## (3 observations deleted due to missingness)
## Multiple R-squared: 0.9374, Adjusted R-squared: 0.9366
## F-statistic: 1184 on 1 and 79 DF, p-value: < 2.2e-16Scatterplot the residuals against the fitted values to check homogeneity of variance
plot(resid(bs2506_FL_lm1) ~ fitted(bs2506_FL_lm1))
The scatterplot clearly looks like it indicates a nonlinear relationship
Plot a histogram to check normality of residuals
qqnorm(bs2506_FL_lm1$residuals)
qqline(bs2506_FL_lm1$residuals)
Calculate AIC
AIC(bs2506_FL_lm1)## [1] 300.5242Plot the model
bs2506_FL_lm1_plottemp <- ggplot(bsdf_2506_FL_W, aes(x = FL_CM, y = W_G)) +
geom_point(aes(FL_CM, W_G),
color = "#242124",
alpha = 0.8) +
labs(x = "Fork Length in cm",
y = "Weight in grams",
title = "Blackfin Snapper (n = 81)") +
theme_classic() +
geom_smooth(method = lm,
aes(FL_CM, W_G),
linewidth = 0.5,
se = FALSE,
color = safe_colorblind_palette[3])
# add label with equation
bs2506_FL_lm1_plot <- bs2506_FL_lm1_plottemp +
geom_text(x = 20, y = 1400,
label = expression(W~'(g)' == 0.018 * FL~'(cm)' + 19.93~"," ~ R^{2} == 0.937 ~ "(Linear)"),
colour = safe_colorblind_palette[3],
hjust = 0)
bs2506_FL_lm1_plot
NLS Models
NLS 1:
bs2506_FL_nls1 <- nls_multstart(W_G~A*FL_CM^B,
data = bsdf_2506_FL_W,
iter = 500,
start_lower = c(A=0, B=2),
start_upper = c(A=0.5, B=3.5))
bs2506_FL_nls1## Nonlinear regression model
## model: W_G ~ A * FL_CM^B
## data: data
## A B
## 0.01517 3.04397
## residual sum-of-squares: 165664
##
## Number of iterations to convergence: 9
## Achieved convergence tolerance: 1.49e-08summary(bs2506_FL_nls1)##
## Formula: W_G ~ A * FL_CM^B
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## A 0.015169 0.002759 5.498 4.57e-07 ***
## B 3.043974 0.050812 59.906 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 45.79 on 79 degrees of freedom
##
## Number of iterations to convergence: 9
## Achieved convergence tolerance: 1.49e-08Plot a histogram of the residuals to check normality of residuals
qqnorm(resid(bs2506_FL_nls1))
qqline(resid(bs2506_FL_nls1))
positive skew
Calculate R-squared
rss_fitted <- sum(residuals(bs2506_FL_nls1)^2)
rss_constant <- sum((bsdf_2506_FL_W$W_G - mean(bsdf_2506_FL_W$W_G))^2)
approx_r_squared <- 1 - rss_fitted / rss_constant
approx_r_squared## [1] NACalculate AIC
AIC(bs2506_FL_nls1)## [1] 853.3528Calculate the predicted W_kg values based on this NLS model
bs2506_FL_nls1_x <- c(18:43) #rounded actual min and max FL_MM values
bs2506_FL_nls1_predic <- data.frame(bs2506_FL_nls1_x)
bs2506_FL_nls1_f <- function(bs2506_FL_nls1_x) {
#formula using calculated parameters
return(0.015169*bs2506_FL_nls1_x^3.043974)
}
bs2506_FL_nls1_predic$bs2506_FL_nls1_y <- sapply(bs2506_FL_nls1_predic$bs2506_FL_nls1_x, bs2506_FL_nls1_f)
View(bs2506_FL_nls1_predic)Plot the model on the same graph as the other model
bs2506_FL_nls1_plottemp <- bs2506_FL_lm1_plot +
geom_smooth(aes(x = bs2506_FL_nls1_x, y = bs2506_FL_nls1_y),
data = bs2506_FL_nls1_predic,
method = "loess",
se = F,
color = safe_colorblind_palette[7],
linewidth = 0.75)
# add label with equation
bs2506_FL_nls1_plot <- bs2506_FL_nls1_plottemp +
geom_text(x = 20, y = 1300,
label = expression(W~'(g)' == 0.015 * TL~'(cm)'^{3.044}~"," ~ R^{2} == 0.980 ~ '(NLS)'),
colour = safe_colorblind_palette[7],
hjust = 0)
bs2506_FL_nls1_plot
Since there is a positive skew in the residuals plot, I will try a few transformations to see if I can find a better fit.
Data transformations
bsdf_2506_FL_W$logFL_CM <- log(bsdf_2506_FL_W$FL_CM)
bsdf_2506_FL_W$logW_G <- log(bsdf_2506_FL_W$W_G)
bsdf_2506_FL_W$log10FL_CM <- log(bsdf_2506_FL_W$FL_CM, 10)
bsdf_2506_FL_W$log10W_G <- log(bsdf_2506_FL_W$W_G, 10)
summary(bsdf_2506_FL_W)## X catch_date catch_month landing_time
## Min. : 2.00 Length:84 : 3 Length:84
## 1st Qu.:22.75 Class :character Apr:10 Class :character
## Median :43.50 Mode :character Dec: 1 Mode :character
## Mean :43.50 Jul: 9
## 3rd Qu.:64.25 Jun:32
## Max. :85.00 May:19
## Sep:10
## dissection_date dissection_time year month
## Length:84 Length:84 2023:52 Length:84
## Class :character Class :character 2024:32 Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
## observer species common_name
## Length:84 Lutjanus buccanella:84 Blackfin Snapper:84
## Class :character
## Mode :character
##
##
##
##
## fish_code frozen catch_method bait
## Length:84 N:32 :51 :51
## Class :character Y:52 Lobster trap:10 Cowskin :10
## Mode :character Redfish trap:23 Japanese bait :22
## Japanese bait; Squid: 1
##
##
##
## area SL_CM FL_CM TL_CM HH_CM
## :41 Min. :16.40 Min. :18.70 Min. :20.50 Min. :3.600
## A3/A4 :17 1st Qu.:24.60 1st Qu.:23.60 1st Qu.:26.00 1st Qu.:6.200
## C2 :10 Median :28.00 Median :26.25 Median :28.30 Median :7.950
## D3 :10 Mean :28.44 Mean :27.67 Mean :29.89 Mean :7.359
## D5 : 3 3rd Qu.:34.20 3rd Qu.:30.05 3rd Qu.:32.50 3rd Qu.:8.850
## D4 : 2 Max. :38.30 Max. :42.90 Max. :45.90 Max. :9.900
## (Other): 1 NA's :51 NA's :3 NA's :62
## ED_CM W_G finclip otoliths gutted
## Min. :1.200 Min. : 108.5 Length:84 Length:84 N:84
## 1st Qu.:2.200 1st Qu.: 234.3 Class :character Class :character Y: 0
## Median :2.400 Median : 323.4 Mode :character Mode :character
## Mean :2.245 Mean : 430.3
## 3rd Qu.:2.500 3rd Qu.: 474.0
## Max. :2.900 Max. :1493.9
## NA's :62 NA's :3
## gonads_present sex maturity G1L_CM G2L_CM
## Length:84 :32 :52 Min. :1.100 Min. :1.200
## Class :character F:18 A:15 1st Qu.:2.500 1st Qu.:2.550
## Mode :character M:34 B:10 Median :3.200 Median :3.200
## C: 7 Mean :3.709 Mean :3.705
## 3rd Qu.:4.200 3rd Qu.:4.500
## Max. :9.300 Max. :8.800
## NA's :19 NA's :25
## GW_G stomach._sample logFL_CM logW_G
## Length:84 Length:84 Min. :2.929 Min. :4.687
## Class :character Class :character 1st Qu.:3.161 1st Qu.:5.457
## Mode :character Mode :character Median :3.268 Median :5.779
## Mean :3.300 Mean :5.856
## 3rd Qu.:3.403 3rd Qu.:6.161
## Max. :3.759 Max. :7.309
## NA's :3
## log10FL_CM log10W_G
## Min. :1.272 Min. :2.035
## 1st Qu.:1.373 1st Qu.:2.370
## Median :1.419 Median :2.510
## Mean :1.433 Mean :2.543
## 3rd Qu.:1.478 3rd Qu.:2.676
## Max. :1.632 Max. :3.174
## NA's :3View(bsdf_2506_FL_W)Natural log-log Linear model (lnFL-lnW)
bs2506_logFL_lm1 <- lm(logFL_CM ~ logW_G, data = bsdf_2506_FL_W)
bs2506_logFL_lm1##
## Call:
## lm(formula = logFL_CM ~ logW_G, data = bsdf_2506_FL_W)
##
## Coefficients:
## (Intercept) logW_G
## 1.4036 0.3237summary(bs2506_logFL_lm1)##
## Call:
## lm(formula = logFL_CM ~ logW_G, data = bsdf_2506_FL_W)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.112947 -0.018656 -0.002307 0.015210 0.163094
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.403634 0.035730 39.28 <2e-16 ***
## logW_G 0.323709 0.006069 53.34 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03346 on 79 degrees of freedom
## (3 observations deleted due to missingness)
## Multiple R-squared: 0.973, Adjusted R-squared: 0.9726
## F-statistic: 2845 on 1 and 79 DF, p-value: < 2.2e-16Scatterplot the residuals against the fitted values to check homogeneity of variance
plot(resid(bs2506_logFL_lm1) ~ fitted(bs2506_logFL_lm1))
Pretty good
Plot a histogram to check normality of residuals
qqnorm(bs2506_logFL_lm1$residuals)
qqline(bs2506_logFL_lm1$residuals)
Okay
Calculate AIC
AIC(bs2506_logFL_lm1)## [1] -316.5224Plot the model
bs2506_logFL_lm1_plottemp <- ggplot(bsdf_2506_FL_W, aes(x = logFL_CM, y = logW_G)) +
geom_point(aes(logFL_CM, logW_G),
color = "#242124",
alpha = 0.8) +
labs(x = "ln(Fork Length in cm)",
y = "ln(Weight in grams)",
title = "Natural Log Transformation of Blackfin Snapper FL-W Relationship (n = 81)") +
theme_classic() +
geom_smooth(method = lm,
aes(logFL_CM, logW_G),
linewidth = 0.5,
se = FALSE,
color = safe_colorblind_palette[3])
# add label with equation
bs2506_logFL_lm1_plot <- bs2506_logFL_lm1_plottemp +
geom_text(x = 20, y = 1400,
label = expression(W~'(g)' == 0.324 * FL~'(cm)' + 1.404~"," ~ R^{2} == 0.937 ~ "(Linear)"),
colour = safe_colorblind_palette[3],
hjust = 0)
bs2506_logFL_lm1_plot
ln NLS:
bs2506_logFL_nls1 <- nls_multstart(logW_G~A*logFL_CM^B,
data = bsdf_2506_FL_W,
iter = 500,
start_lower = c(A=0, B=1),
start_upper = c(A=1, B=3))
bs2506_logFL_nls1## Nonlinear regression model
## model: logW_G ~ A * logFL_CM^B
## data: data
## A B
## 0.7934 1.6728
## residual sum-of-squares: 0.8233
##
## Number of iterations to convergence: 6
## Achieved convergence tolerance: 1.49e-08summary(bs2506_logFL_nls1)##
## Formula: logW_G ~ A * logFL_CM^B
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## A 0.79341 0.02934 27.05 <2e-16 ***
## B 1.67279 0.03066 54.57 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1021 on 79 degrees of freedom
##
## Number of iterations to convergence: 6
## Achieved convergence tolerance: 1.49e-08A = 0.79341 B = 1.67279
Plot a histogram of the residuals to check normality of residuals
qqnorm(resid(bs2506_logFL_nls1))
qqline(resid(bs2506_logFL_nls1))
Calculate R-squared
rss_fitted <- sum(residuals(bs2506_logFL_nls1)^2)
rss_constant <- sum((bsdf_2506_FL_W$W_G - mean(bsdf_2506_FL_W$W_G))^2)
approx_r_squared <- 1 - rss_fitted / rss_constant
approx_r_squared## [1] NACalculate AIC
AIC(bs2506_logFL_nls1)## [1] -135.8284Calculate the predicted W_kg values based on this NLS model
bs2506_logFL_nls1_x <- seq(2.9,3.8, by = 0.1) #rounded actual min and max FL_MM values
bs2506_logFL_nls1_predic <- data.frame(bs2506_logFL_nls1_x)
bs2506_logFL_nls1_f <- function(bs2506_logFL_nls1_x) {
#formula using calculated parameters
return(0.79341*bs2506_logFL_nls1_x^1.67279)
}
bs2506_logFL_nls1_predic$bs2506_logFL_nls1_y <- sapply(bs2506_logFL_nls1_predic$bs2506_logFL_nls1_x, bs2506_logFL_nls1_f)
View(bs2506_logFL_nls1_predic)Plot the model on the same graph as the other model
bs2506_logFL_nls1_plottemp <- bs2506_logFL_lm1_plot +
geom_smooth(aes(x = bs2506_logFL_nls1_x, y = bs2506_logFL_nls1_y),
data = bs2506_logFL_nls1_predic,
method = "loess",
se = F,
color = safe_colorblind_palette[7],
linewidth = 0.75,
linetype = 2)
# add label with equation
bs2506_logFL_nls1_plot <- bs2506_logFL_nls1_plottemp +
geom_text(x = 20, y = 1300,
label = expression(W~'(g)' == 0.793 * TL~'(cm)'^{1.673}~"," ~ R^{2} == 0.999 ~ '(NLS)'),
colour = safe_colorblind_palette[7],
hjust = 0)
bs2506_logFL_nls1_plot## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 3 rows containing non-finite outside the scale range
## (`stat_smooth()`).## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'
## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'
Log10 Linear model *
bs2506_log10FL_lm1 <- lm(log10FL_CM ~ log10W_G, data = bsdf_2506_FL_W)
summary(bs2506_log10FL_lm1)##
## Call:
## lm(formula = log10FL_CM ~ log10W_G, data = bsdf_2506_FL_W)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.049052 -0.008102 -0.001002 0.006605 0.070831
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.609591 0.015518 39.28 <2e-16 ***
## log10W_G 0.323709 0.006069 53.34 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01453 on 79 degrees of freedom
## (3 observations deleted due to missingness)
## Multiple R-squared: 0.973, Adjusted R-squared: 0.9726
## F-statistic: 2845 on 1 and 79 DF, p-value: < 2.2e-16length(bs2506_log10FL_lm1$fitted.values)## [1] 81Scatterplot the residuals against the fitted values to check homogeneity of variance
plot(resid(bs2506_log10FL_lm1) ~ fitted(bs2506_log10FL_lm1))
Pretty good
Plot a histogram to check normality of residuals
qqnorm(bs2506_log10FL_lm1$residuals)
qqline(bs2506_log10FL_lm1$residuals)
fine
Calculate AIC
AIC(bs2506_log10FL_lm1)## [1] -451.6357Plot the model
bs2506_log10FL_lm1_plot_title <- expression(Log[10]-Log[10]~"Transformation of"~italic("L. buccanella")~"FL-W Relationship (n = 81)")
bs2506_log10FL_lm1_plottemp <- ggplot(bsdf_2506_FL_W, aes(x = log10FL_CM, y = log10W_G)) +
geom_point(aes(log10FL_CM, log10W_G),
color = "#242124",
alpha = 0.8) +
labs(x = expression(log[10]~"Fork Length in cm"),
y = expression(log[10]~"Weight in grams"),
title = bs2506_log10FL_lm1_plot_title) +
theme_classic() +
geom_smooth(method = lm,
aes(log10FL_CM, log10W_G),
linewidth = 0.75,
se = FALSE,
color = turbo_colors_onboard[3]) +
theme(plot.margin = unit(c(0.5, 1, 0.5, 0.5), "cm"))
# add label with equation
bs2506_log10FL_lm1_plot <- bs2506_log10FL_lm1_plottemp +
geom_label(x = 1.27, y = 3.3,
label = expression(log[10] ~ "W(g)" == 0.324 * Log[10] ~ "FL(cm)" + 0.610~"," ~ r^{2} == 0.973~"(Linear)"),
fill = "#91e4f0",
hjust = 0)
bs2506_log10FL_lm1_plot
NLS Models
NLS 1:
bs2506_log10FL_nls1 <- nls_multstart(log10W_G~A*log10FL_CM^B,
data = bsdf_2506_FL_W,
iter = 500,
start_lower = c(A=0, B=1),
start_upper = c(A=1, B=3))
bs2506_log10FL_nls1## Nonlinear regression model
## model: log10W_G ~ A * log10FL_CM^B
## data: data
## A B
## 1.391 1.673
## residual sum-of-squares: 0.1553
##
## Number of iterations to convergence: 7
## Achieved convergence tolerance: 1.49e-08summary(bs2506_log10FL_nls1)##
## Formula: log10W_G ~ A * log10FL_CM^B
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## A 1.39056 0.01602 86.81 <2e-16 ***
## B 1.67279 0.03066 54.57 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04434 on 79 degrees of freedom
##
## Number of iterations to convergence: 7
## Achieved convergence tolerance: 1.49e-08Plot a histogram of the residuals to check normality of residuals
qqnorm(resid(bs2506_log10FL_nls1))
qqline(resid(bs2506_log10FL_nls1))
Calculate R-squared
rss_fitted <- sum(residuals(bs2506_log10FL_nls1)^2)
rss_constant <- sum((bsdf_2506_FL_W$W_G - mean(bsdf_2506_FL_W$W_G))^2)
approx_r_squared <- 1 - rss_fitted / rss_constant
approx_r_squared## [1] NACalculate AIC
AIC(bs2506_log10FL_nls1)## [1] -270.9417Calculate the predicted W_kg values based on this NLS model
bs2506_log10FL_nls1_x <- seq(1.27,1.64, by = 0.01) #rounded actual min and max FL_MM values
bs2506_log10FL_nls1_predic <- data.frame(bs2506_log10FL_nls1_x)
bs2506_log10FL_nls1_f <- function(bs2506_log10FL_nls1_x) {
#formula using calculated parameters
return(1.39056*bs2506_log10FL_nls1_x^1.67279)
}
bs2506_log10FL_nls1_predic$bs2506_log10FL_nls1_y <- sapply(bs2506_log10FL_nls1_predic$bs2506_log10FL_nls1_x, bs2506_log10FL_nls1_f)
View(bs2506_log10FL_nls1_predic)Plot the model on the same graph as the other model
bs2506_log10FL_nls1_plottemp <- bs2506_log10FL_lm1_plot +
geom_smooth(aes(x = bs2506_log10FL_nls1_x, y = bs2506_log10FL_nls1_y),
data = bs2506_log10FL_nls1_predic,
method = "loess",
se = F,
color = safe_colorblind_palette[7],
linewidth = 0.75,
linetype = 2)
# add label with equation
bs2506_log10FL_nls1_plot <- bs2506_log10FL_nls1_plottemp +
geom_text(x = 1.27, y = 3,
label = expression(W~'(g)' == 1.391 * TL~'(cm)'^{1.673}~"," ~ R^{2} == 1 ~ '(NLS)'),
colour = safe_colorblind_palette[7],
hjust = 0)
bs2506_log10FL_nls1_plot## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 3 rows containing non-finite outside the scale range
## (`stat_smooth()`).## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'
## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'
#Final decision: Log10-Log10 Linear Model
Because there is a negative skew in the residuals, I think the linear model fits this relationship better.
Save plot
ggsave("bs2506_FLlog10_lm1_plot.png", plot = bs2506_log10FL_lm1_plot, width = 6, height = 4, dpi = 400)## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 3 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'Total length (TL_CM) - Weight (W_G) relationship
Load data file
bsdf_2506_TL_W <- read.csv("bs2506df.csv")
bsdf_2506_TL_W$catch_month <- factor(bsdf_2506_TL_W$catch_month)
bsdf_2506_TL_W$year <- factor(bsdf_2506_TL_W$year)
bsdf_2506_TL_W$species <- factor(bsdf_2506_TL_W$species)
bsdf_2506_TL_W$common_name <- factor(bsdf_2506_TL_W$common_name)
bsdf_2506_TL_W$frozen <- factor(bsdf_2506_TL_W$frozen)
bsdf_2506_TL_W$catch_method <- factor(bsdf_2506_TL_W$catch_method)
bsdf_2506_TL_W$bait <- factor(bsdf_2506_TL_W$bait)
bsdf_2506_TL_W$area <- factor(bsdf_2506_TL_W$area)
bsdf_2506_TL_W$gutted <- factor(bsdf_2506_TL_W$gutted)
bsdf_2506_TL_W$sex <- factor(bsdf_2506_TL_W$sex)
bsdf_2506_TL_W$maturity <- factor(bsdf_2506_TL_W$maturity)
# Remove gutted samples
bsdf_2506_TL_W <- subset(filter(bsdf_2506_TL_W, gutted == "N"))
View(bsdf_2506_TL_W)
summary(bsdf_2506_TL_W)## X catch_date catch_month landing_time
## Min. : 2.00 Length:84 : 3 Length:84
## 1st Qu.:22.75 Class :character Apr:10 Class :character
## Median :43.50 Mode :character Dec: 1 Mode :character
## Mean :43.50 Jul: 9
## 3rd Qu.:64.25 Jun:32
## Max. :85.00 May:19
## Sep:10
## dissection_date dissection_time year month
## Length:84 Length:84 2023:52 Length:84
## Class :character Class :character 2024:32 Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
## observer species common_name
## Length:84 Lutjanus buccanella:84 Blackfin Snapper:84
## Class :character
## Mode :character
##
##
##
##
## fish_code frozen catch_method bait
## Length:84 N:32 :51 :51
## Class :character Y:52 Lobster trap:10 Cowskin :10
## Mode :character Redfish trap:23 Japanese bait :22
## Japanese bait; Squid: 1
##
##
##
## area SL_CM FL_CM TL_CM HH_CM
## :41 Min. :16.40 Min. :18.70 Min. :20.50 Min. :3.600
## A3/A4 :17 1st Qu.:24.60 1st Qu.:23.60 1st Qu.:26.00 1st Qu.:6.200
## C2 :10 Median :28.00 Median :26.25 Median :28.30 Median :7.950
## D3 :10 Mean :28.44 Mean :27.67 Mean :29.89 Mean :7.359
## D5 : 3 3rd Qu.:34.20 3rd Qu.:30.05 3rd Qu.:32.50 3rd Qu.:8.850
## D4 : 2 Max. :38.30 Max. :42.90 Max. :45.90 Max. :9.900
## (Other): 1 NA's :51 NA's :3 NA's :62
## ED_CM W_G finclip otoliths gutted
## Min. :1.200 Min. : 108.5 Length:84 Length:84 N:84
## 1st Qu.:2.200 1st Qu.: 234.3 Class :character Class :character Y: 0
## Median :2.400 Median : 323.4 Mode :character Mode :character
## Mean :2.245 Mean : 430.3
## 3rd Qu.:2.500 3rd Qu.: 474.0
## Max. :2.900 Max. :1493.9
## NA's :62 NA's :3
## gonads_present sex maturity G1L_CM G2L_CM
## Length:84 :32 :52 Min. :1.100 Min. :1.200
## Class :character F:18 A:15 1st Qu.:2.500 1st Qu.:2.550
## Mode :character M:34 B:10 Median :3.200 Median :3.200
## C: 7 Mean :3.709 Mean :3.705
## 3rd Qu.:4.200 3rd Qu.:4.500
## Max. :9.300 Max. :8.800
## NA's :19 NA's :25
## GW_G stomach._sample
## Length:84 Length:84
## Class :character Class :character
## Mode :character Mode :character
##
##
##
## First look at the data with a scatterplot
fill = sex
ggplot(bsdf_2506_TL_W, aes(x = TL_CM, y = W_G)) +
geom_point(aes(fill = sex), color = "#1B1B1B", pch = 21, size = 2) +
theme_classic() +
labs(x = "Total Length in cm", y = "Weight in grams", title = "Lutjanus buccanella, n = 80") +
scale_fill_viridis(option = "D", discrete = TRUE)
Curves from published literature on L. buccanella (BS) TL-W relationships
W_kg values calculated from Burton et al. (2016) parameters then converted to cm and g
Use size range from their paper (TL = 180mm-609mm) W = 7.79e-09TL^(3.09)
#predicted values
bs_burton_x <- c(180:609)
bs_burton_df <- data.frame(bs_burton_x)
bs_burton_f <- function(bs_burton_x) {
# Replace this with your actual formula
return(7.79e-9*bs_burton_x^3.09)
}
bs_burton_df$bs_burton_y <- sapply(bs_burton_df$bs_burton_x, bs_burton_f)
#change kg to g
bs_burton_df$bs_burton_y_g <- bs_burton_df$bs_burton_y*1000
#change mm to cm
bs_burton_df$bs_burton_x_cm <- bs_burton_df$bs_burton_x/10
View(bs_burton_df)W_g values calculated from Tabash and Sierra (1996) parameters
On natural log transformed data W(g) = 0.0000142 L(mm) ^ 2.89 (N = 200, rsq = 0.864) ###### CAN’T COMPARE ON NORMAL GRAPH - make ln transformed graph
bs_tabash_x <- seq(from = log(180), to = log(540), by = 0.01)
bs_tabash_df <- data.frame(bs_tabash_x)
bs_tabash_f <- function(bs_tabash_x) {
return(0.0000142*bs_tabash_x^2.89)
}
bs_tabash_df$y <- sapply(bs_tabash_df$bs_tabash_x, bs_tabash_f)
#change log(TL) to TL_mm
bs_tabash_df$bs_tabash_x_mm <- exp(bs_tabash_df$bs_tabash_x)
#change mm to cm
bs_tabash_df$bs_tabash_x_cm <- bs_tabash_df$bs_tabash_x_mm/10
#change log(W) to W_g
bs_tabash_df$y_g <- exp(bs_tabash_df$y)
View(bs_tabash_df)Plot models from literature
bs_lit_plottemp <- ggplot(bsdf_2506_TL_W, aes(x = TL_CM, y = W_G)) +
geom_point(aes(TL_CM, W_G),
color = "#242124",
alpha = 0.8) +
geom_smooth(aes(x = bs_burton_x_cm, y = bs_burton_y_g),
data = bs_burton_df,
method = "loess",
se = FALSE,
linewidth = 0.75,
linetype = 2,
color = safe_colorblind_palette[1]) +
labs(x = "Total Length, cm",
y = "Weight, grams",
title = "Blackfin Snapper: Saba vs. Other Regions (n = 84)") +
theme_classic()
# add label with equation 7.79e-9*bs_burton_x^3.09
bs_lit_plot <- bs_lit_plottemp +
geom_text(x = 18, y = 2800,
label = expression(W == 7.79 %*% 10^{-09} * TL^{3.09} ~ "(SE & Caribbean U.S.)"),
colour = safe_colorblind_palette[1],
hjust = 0)
bs_lit_plot
Linear model
remove NA rows
bsdf_2506_TL_W <- bsdf_2506_TL_W%>%
drop_na(TL_CM)
View(bsdf_2506_TL_W)bs2506_lm1 <- lm(TL_CM ~ W_G, data = bsdf_2506_TL_W)
summary(bs2506_lm1)##
## Call:
## lm(formula = TL_CM ~ W_G, data = bsdf_2506_TL_W)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.8780 -0.9322 0.2306 1.2948 2.9979
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.173e+01 3.014e-01 72.07 <2e-16 ***
## W_G 1.898e-02 5.616e-04 33.79 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.622 on 79 degrees of freedom
## Multiple R-squared: 0.9353, Adjusted R-squared: 0.9345
## F-statistic: 1142 on 1 and 79 DF, p-value: < 2.2e-16Scatterplot the residuals against the fitted values to check homogeneity of variance
plot(resid(bs2506_lm1) ~ fitted(bs2506_lm1))
The scatterplot clearly looks like it indicates a nonlinear relationship
Plot a histogram to check normality of residuals
qqnorm(resid(bs2506_lm1))
qqline(resid(bs2506_lm1))
negative skew
Calculate AIC
AIC(bs2506_lm1)## [1] 312.1639Plot the model
bs2506_lm1_plottemp <- ggplot(bsdf_2506_TL_W, aes(x = TL_CM, y = W_G)) +
geom_point(aes(TL_CM, W_G),
color = "#242124",
alpha = 0.8) +
labs(x = "Total Length in cm",
y = "Weight in grams",
title = "Blackfin Snapper (n = 81)") +
theme_classic() +
geom_smooth(method = lm,
aes(TL_CM, W_G),
linewidth = 0.5,
se = FALSE,
color = safe_colorblind_palette[3])
# add label with equation
bs2506_lm1_plot <- bs2506_lm1_plottemp +
geom_text(x = 20, y = 1400,
label = expression(W~'(g)' == 0.020 * TL~'(cm)' + 21.479~"," ~ R^{2} == 0.935 ~ "(Linear)"),
colour = safe_colorblind_palette[3],
hjust = 0)
bs2506_lm1_plot
NLS Models
NLS 1:
bs2506_nls1 <- nls_multstart(W_G~A*TL_CM^B,
data = bsdf_2506_TL_W,
iter = 500,
start_lower = c(A=0, B=2),
start_upper = c(A=0.5, B=3.5))
summary(bs2506_nls1)##
## Formula: W_G ~ A * TL_CM^B
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## A 0.009687 0.001869 5.184 1.63e-06 ***
## B 3.106410 0.052825 58.806 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 46.7 on 79 degrees of freedom
##
## Number of iterations to convergence: 13
## Achieved convergence tolerance: 1.49e-08Plot a histogram of the residuals to check normality of residuals
qqnorm(resid(bs2506_nls1))
qqline(resid(bs2506_nls1))
heavy tailed but okay
Calculate R-squared
rss_fitted <- sum(residuals(bs2506_nls1)^2)
rss_constant <- sum((bsdf_2506_TL_W$W_G - mean(bsdf_2506_TL_W$W_G))^2)
approx_r_squared <- 1 - rss_fitted / rss_constant
approx_r_squared## [1] 0.9793393Calculate AIC
AIC(bs2506_nls1)## [1] 856.5134Calculate the predicted W_kg values based on this NLS model
bs2506_nls1_x <- c(20:46) #rounded actual min and max TL_CM values
bs2506_nls1_predic <- data.frame(bs2506_nls1_x)
bs2506_nls1_f <- function(bs2506_nls1_x) {
#formula using calculated parameters
return(0.009687*bs2506_nls1_x^3.106410)
}
bs2506_nls1_predic$bs2506_nls1_y <- sapply(bs2506_nls1_predic$bs2506_nls1_x, bs2506_nls1_f)
View(bs2506_nls1_predic)Plot the model on the same graph as the other model
bs2506_nls1_plottemp <- bs2506_lm1_plot +
geom_smooth(aes(x = bs2506_nls1_x, y = bs2506_nls1_y),
data = bs2506_nls1_predic,
method = "loess",
se = F,
color = safe_colorblind_palette[7],
linewidth = 0.75)
# add label with equation
bs2506_nls1_plot <- bs2506_nls1_plottemp +
geom_text(x = 20, y = 1300,
label = expression(W~'(g)' == 9.687%*%10^{-03} * TL~'(cm)'^{3.106}~"," ~ R^{2} == 0.979 ~ '(NLS)'),
colour = safe_colorblind_palette[7],
hjust = 0)
bs2506_nls1_plot
NLS2: model starting values based on normal distribution around parameters from previous L. buccanella paper (Burton et al., 2016)
W = A*TL^B, n = 216 A = 7.79e-9) (SE = 4.47e-9) B = 3.09 (SE= 0.09)
SD_A = 6.569531e-08 SD_B = 1.322724
Need to change the units of measurement to compare to Burton et al. (2016) which used mm and kg
TL_MM <- bsdf_2506_TL_W$TL_CM*10
bsdf_2506_TL_W$TL_MM <- TL_MM
W_kg <- bsdf_2506_TL_W$W_G/1000
bsdf_2506_TL_W$W_kg <- W_kg
View(bsdf_2506_TL_W)Run the model
bs2506_nls2 <- nls_multstart(W_kg~A*TL_MM^B,
data = bsdf_2506_TL_W,
iter = 500,
start_lower = c(A=0, B=1.232724),
start_upper = c(A=0.001, B=3.18)) # upper and lower bounds for B are calculated with the SE of the mean
summary(bs2506_nls2)##
## Formula: W_kg ~ A * TL_MM^B
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## A 7.582e-09 2.384e-09 3.18 0.0021 **
## B 3.106e+00 5.283e-02 58.81 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0467 on 79 degrees of freedom
##
## Number of iterations to convergence: 60
## Achieved convergence tolerance: 1.49e-08The calculated parameter for A is smaller than the previous NLS model but the B value is the same
Plot a histogram of the residuals to check normality of residuals
qqnorm(resid(bs2506_nls2))
qqline(resid(bs2506_nls2))
Calculate R-squared
rss_fitted <- sum(residuals(bs2506_nls2)^2)
rss_constant <- sum((bsdf_2506_TL_W$W_kg - mean(bsdf_2506_TL_W$W_kg))^2)
approx_r_squared <- 1 - rss_fitted / rss_constant
approx_r_squared ## [1] 0.9793393Calculate AIC
AIC(bs2506_nls2)## [1] -262.5429Calculate the predicted W_kg values based on this NLS model
bs2506_nls2_x <- c(200:460) #rounded actual min and max TL_MM values
bs2506_nls2_predic <- data.frame(bs2506_nls2_x)
bs2506_nls2_f <- function(bs2506_nls2_x) {
#formula using calculated parameters
return(7.582e-09*bs2506_nls2_x^3.106) }
bs2506_nls2_predic$bs2506_nls2_y <- sapply(bs2506_nls2_predic$bs2506_nls2_x, bs2506_nls2_f)
View(bs2506_nls2_predic)Convert predicted values back from TL_MM to TL_CM and W_kg to W_G
bs2506_nls2_predic$bs2506_nls2_x_CM <- bs2506_nls2_x/10
bs2506_nls2_predic$bs2506_nls2_y_g <- bs2506_nls2_predic$bs2506_nls2_y*1000
View(bs2506_nls2_predic)Plot the model on the same graph as the other model
bs2506_nls2_plottemp <- bs2506_nls1_plot +
geom_smooth(aes(x = bs2506_nls2_x_CM, y = bs2506_nls2_y_g),
data = bs2506_nls2_predic,
method = "loess",
se = F,
color = safe_colorblind_palette[6],
linewidth = 1,
linetype = 2)
# add label with equation
bs2506_nls2_plot <- bs2506_nls2_plottemp +
geom_text(x = 20, y = 1200,
label = expression(W~'(kg)' == 7.582%*%10^{-09} * TL~'(mm)'^{3.106}~"," ~ R^{2} == 0.979 ~ '(NLS)'),
colour = safe_colorblind_palette[6],
hjust = 0)
bs2506_nls2_plot
Since these parameter values seem to fit the data well according to the r-squared and the AIC value, I will use them as a starting point for Bayesian modeling.
Bayesian models
Brms1
Very small prior values often do not work, so for A I will choose a slightly larger range than the exact value predicted in NLS
bs2506_prior1 <- prior(normal(0.01, 0.5), nlpar = "A") +
prior(normal(3.106, 1), nlpar = "B")
#brm doesn't seem to work if the conflicted package is installed so restart R if it is installed and don't load it
bs2506_brm1 <- brm(bf(W_kg ~ A * TL_MM ^ B, A + B ~ 1, nl = TRUE),
data = bsdf_2506_TL_W, prior = bs2506_prior1,
warmup = 1000,
iter = 6000,
chains = 4,
cores = 4,
control = list(adapt_delta = 0.8, max_treedepth = 12))## Compiling Stan program...## Start sampling## Warning: There were 786 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.## Warning: There were 17282 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 12. See
## https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded## Warning: Examine the pairs() plot to diagnose sampling problems## Warning: The largest R-hat is 1.84, indicating chains have not mixed.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#r-hat## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-essbs2506_brm1## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.## Warning: There were 786 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: W_kg ~ A * TL_MM^B
## A ~ 1
## B ~ 1
## Data: bsdf_2506_TL_W (Number of observations: 81)
## Draws: 4 chains, each with iter = 6000; warmup = 1000; thin = 1;
## total post-warmup draws = 20000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## A_Intercept 0.00 0.00 0.00 0.00 1.83 6 16
## B_Intercept 2.15 0.06 2.06 2.26 1.84 6 22
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.15 0.01 0.12 0.16 1.75 6 15
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).Lots of warnings on this model and A = 0 which will result in all W values equaling 0. I’ll try going back a step and use the actual Burton values as priors for B and a wider estimate for A.
BRM2
Adjust so mu is bigger than prior value
bs2506_prior2 <- prior(normal(0.03, 0.01), nlpar = "A") +
prior(normal(3.106, 1), nlpar = "B")
#brm doesn't seem to work if the conflicted package is installed so restart R if it is installed and don't load it
bs2506_brm2 <- brm(bf(W_G ~ A * TL_CM ^ B, A + B ~ 1, nl = TRUE),
data = bsdf_2506_TL_W, prior = bs2506_prior2,
warmup = 1000,
iter = 6000,
chains = 4,
cores = 4,
control = list(adapt_delta = 0.8))## Compiling Stan program...## Start samplingprint(summary(bs2506_brm2), digits = 6)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: W_G ~ A * TL_CM^B
## A ~ 1
## B ~ 1
## Data: bsdf_2506_TL_W (Number of observations: 81)
## Draws: 4 chains, each with iter = 6000; warmup = 1000; thin = 1;
## total post-warmup draws = 20000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## A_Intercept 0.011089 0.002191 0.007427 0.015970 1.000976 5425 5901
## B_Intercept 3.074554 0.053338 2.968558 3.179087 1.000810 5429 5938
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 47.453410 3.813459 40.636442 55.525964 1.000185 6931 7795
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).Calculate R-squared
rss_fitted <- sum(residuals(bs2506_brm2)^2)
rss_constant <- sum((bsdf_2506_TL_W$W_G - mean(bsdf_2506_TL_W$W_G))^2)
approx_r_squared <- 1 - rss_fitted / rss_constant
approx_r_squared ## [1] 0.7411596bs2506_brm2_x <- c(20:46) #rounded actual min and max TL_MM values
bs2506_brm2_predic <- data.frame(bs2506_brm2_x)
bs2506_brm2_f <- function(bs2506_brm2_x) {
#formula using calculated parameters
return(0.011015*bs2506_brm2_x^3.076398) }
bs2506_brm2_predic$bs2506_brm2_y <- sapply(bs2506_brm2_predic$bs2506_brm2_x, bs2506_brm2_f)
View(bs2506_brm2_predic)bs2506_brm2_plottemp <- bs2506_nls2_plot +
geom_smooth(aes(x = bs2506_brm2_x, y = bs2506_brm2_y),
data = bs2506_brm2_predic,
method = "loess",
se = F,
color = safe_colorblind_palette[8],
linewidth = 1,
linetype = 3)
# add label with equation
bs2506_brm2_plot <- bs2506_brm2_plottemp +
geom_text(x = 20, y = 1100,
label = expression(W~'(g)' == 1.1015%*%10^{-02} * TL~'(cm)'^{3.076}~"," ~ R^{2} == 0.740 ~ '(BRMS)'),
colour = safe_colorblind_palette[8],
hjust = 0)
bs2506_brm2_plot
Final model TL-W:NLS1
Because the brms does not work well for this model, I will use one of the NLS models as the final. As the NLS models are very similar with similar r-squared values, I will use the one that has the same units as my own study (grams and cm) = NLS1.
Plot the model
bs2506_WTLnls1_plot_title <- expression(italic("L. buccanella")~"TL-W Relationship (n = 81)")
bs2506_WTLnls1_plot_temp <- ggplot(bsdf_2506_TL_W, aes(x = TL_CM, y = W_G)) +
geom_point(aes(TL_CM, W_G),
color = "#242124",
alpha = 0.8) +
labs(x = expression("Total Length in cm"),
y = expression("Weight in grams"),
title = bs2506_WTLnls1_plot_title) +
geom_smooth(aes(x = bs2506_nls1_x, y = bs2506_nls1_y),
data = bs2506_nls1_predic,
method = "loess",
se = F,
color = "#1bb2c4",
linewidth = 0.75) +
theme_classic()
bs2506_WTLnls1_plot <- bs2506_WTLnls1_plot_temp +
geom_text(x = 20, y = 1400,
label = expression(W~'(g)' == 9.687%*%10^{-03} * TL~'(cm)'^{3.106}~"," ~ R^{2} == 0.979 ~ '(NLS)'),
colour = "#1bb2c4",
hjust = 0)
bs2506_WTLnls1_plot
look into literature to see if there is anything that says about brms
Standard length (SL_CM) - Weight (W_G) relationship
Load data file
bsdf_2506_SL_W <- read.csv("bs2506df.csv")
bsdf_2506_SL_W$catch_month <- factor(bsdf_2506_SL_W$catch_month)
bsdf_2506_SL_W$year <- factor(bsdf_2506_SL_W$year)
bsdf_2506_SL_W$species <- factor(bsdf_2506_SL_W$species)
bsdf_2506_SL_W$common_name <- factor(bsdf_2506_SL_W$common_name)
bsdf_2506_SL_W$frozen <- factor(bsdf_2506_SL_W$frozen)
bsdf_2506_SL_W$catch_method <- factor(bsdf_2506_SL_W$catch_method)
bsdf_2506_SL_W$bait <- factor(bsdf_2506_SL_W$bait)
bsdf_2506_SL_W$area <- factor(bsdf_2506_SL_W$area)
bsdf_2506_SL_W$gutted <- factor(bsdf_2506_SL_W$gutted)
bsdf_2506_SL_W$sex <- factor(bsdf_2506_SL_W$sex)
bsdf_2506_SL_W$maturity <- factor(bsdf_2506_SL_W$maturity)
# Remove gutted samples
bsdf_2506_SL_W <- subset(filter(bsdf_2506_SL_W, gutted == "N"))
View(bsdf_2506_SL_W)
summary(bsdf_2506_SL_W)## X catch_date catch_month landing_time
## Min. : 2.00 Length:84 : 3 Length:84
## 1st Qu.:22.75 Class :character Apr:10 Class :character
## Median :43.50 Mode :character Dec: 1 Mode :character
## Mean :43.50 Jul: 9
## 3rd Qu.:64.25 Jun:32
## Max. :85.00 May:19
## Sep:10
## dissection_date dissection_time year month
## Length:84 Length:84 2023:52 Length:84
## Class :character Class :character 2024:32 Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
## observer species common_name
## Length:84 Lutjanus buccanella:84 Blackfin Snapper:84
## Class :character
## Mode :character
##
##
##
##
## fish_code frozen catch_method bait
## Length:84 N:32 :51 :51
## Class :character Y:52 Lobster trap:10 Cowskin :10
## Mode :character Redfish trap:23 Japanese bait :22
## Japanese bait; Squid: 1
##
##
##
## area SL_CM FL_CM TL_CM HH_CM
## :41 Min. :16.40 Min. :18.70 Min. :20.50 Min. :3.600
## A3/A4 :17 1st Qu.:24.60 1st Qu.:23.60 1st Qu.:26.00 1st Qu.:6.200
## C2 :10 Median :28.00 Median :26.25 Median :28.30 Median :7.950
## D3 :10 Mean :28.44 Mean :27.67 Mean :29.89 Mean :7.359
## D5 : 3 3rd Qu.:34.20 3rd Qu.:30.05 3rd Qu.:32.50 3rd Qu.:8.850
## D4 : 2 Max. :38.30 Max. :42.90 Max. :45.90 Max. :9.900
## (Other): 1 NA's :51 NA's :3 NA's :62
## ED_CM W_G finclip otoliths gutted
## Min. :1.200 Min. : 108.5 Length:84 Length:84 N:84
## 1st Qu.:2.200 1st Qu.: 234.3 Class :character Class :character Y: 0
## Median :2.400 Median : 323.4 Mode :character Mode :character
## Mean :2.245 Mean : 430.3
## 3rd Qu.:2.500 3rd Qu.: 474.0
## Max. :2.900 Max. :1493.9
## NA's :62 NA's :3
## gonads_present sex maturity G1L_CM G2L_CM
## Length:84 :32 :52 Min. :1.100 Min. :1.200
## Class :character F:18 A:15 1st Qu.:2.500 1st Qu.:2.550
## Mode :character M:34 B:10 Median :3.200 Median :3.200
## C: 7 Mean :3.709 Mean :3.705
## 3rd Qu.:4.200 3rd Qu.:4.500
## Max. :9.300 Max. :8.800
## NA's :19 NA's :25
## GW_G stomach._sample
## Length:84 Length:84
## Class :character Class :character
## Mode :character Mode :character
##
##
##
## # remove NA rows
bsdf_2506_SL_W <- bsdf_2506_SL_W%>%
drop_na(SL_CM)
View(bsdf_2506_SL_W)
summary(bsdf_2506_SL_W)## X catch_date catch_month landing_time
## Min. : 2 Length:33 : 0 Length:33
## 1st Qu.:10 Class :character Apr:10 Class :character
## Median :18 Mode :character Dec: 1 Mode :character
## Mean :18 Jul: 0
## 3rd Qu.:26 Jun: 3
## Max. :34 May:19
## Sep: 0
## dissection_date dissection_time year month
## Length:33 Length:33 2023: 1 Length:33
## Class :character Class :character 2024:32 Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
## observer species common_name
## Length:33 Lutjanus buccanella:33 Blackfin Snapper:33
## Class :character
## Mode :character
##
##
##
##
## fish_code frozen catch_method bait
## Length:33 N:32 : 0 : 0
## Class :character Y: 1 Lobster trap:10 Cowskin :10
## Mode :character Redfish trap:23 Japanese bait :22
## Japanese bait; Squid: 1
##
##
##
## area SL_CM FL_CM TL_CM HH_CM
## A3/A4 :17 Min. :16.40 Min. :18.70 Min. :20.50 Min. :3.600
## D3 :10 1st Qu.:24.60 1st Qu.:27.60 1st Qu.:29.80 1st Qu.:6.200
## D5 : 3 Median :28.00 Median :31.30 Median :33.40 Median :7.950
## D4 : 2 Mean :28.44 Mean :32.28 Mean :34.67 Mean :7.359
## B4 : 1 3rd Qu.:34.20 3rd Qu.:38.90 3rd Qu.:41.50 3rd Qu.:8.850
## : 0 Max. :38.30 Max. :42.90 Max. :45.90 Max. :9.900
## (Other): 0 NA's :11
## ED_CM W_G finclip otoliths gutted
## Min. :1.200 Min. : 108.5 Length:33 Length:33 N:33
## 1st Qu.:2.200 1st Qu.: 352.5 Class :character Class :character Y: 0
## Median :2.400 Median : 492.7 Mode :character Mode :character
## Mean :2.245 Mean : 658.9
## 3rd Qu.:2.500 3rd Qu.:1011.7
## Max. :2.900 Max. :1493.9
## NA's :11
## gonads_present sex maturity G1L_CM G2L_CM
## Length:33 : 2 : 1 Min. :1.600 Min. :1.200
## Class :character F: 7 A:15 1st Qu.:2.950 1st Qu.:2.700
## Mode :character M:24 B:10 Median :3.550 Median :3.500
## C: 7 Mean :3.962 Mean :3.983
## 3rd Qu.:4.500 3rd Qu.:5.050
## Max. :8.500 Max. :8.400
## NA's :1 NA's :3
## GW_G stomach._sample
## Length:33 Length:33
## Class :character Class :character
## Mode :character Mode :character
##
##
##
## Linear model
bs2506_SLW_lm <- lm(SL_CM ~ W_G, data = bsdf_2506_SL_W)
summary(bs2506_SLW_lm)##
## Call:
## lm(formula = SL_CM ~ W_G, data = bsdf_2506_SL_W)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.2163 -0.7104 0.2969 1.1438 2.0210
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.907e+01 5.517e-01 34.57 <2e-16 ***
## W_G 1.422e-02 7.216e-04 19.70 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.608 on 31 degrees of freedom
## Multiple R-squared: 0.9261, Adjusted R-squared: 0.9237
## F-statistic: 388.2 on 1 and 31 DF, p-value: < 2.2e-16Scatterplot the residuals against the fitted values to check homogeneity of variance
plot(resid(bs2506_SLW_lm) ~ fitted(bs2506_SLW_lm))
The scatterplot clearly looks like it indicates a nonlinear relationship
Plot a histogram to check normality of residuals
qqnorm(resid(bs2506_SLW_lm))
qqline(resid(bs2506_SLW_lm))
negative skew
Calculate AIC
AIC(bs2506_SLW_lm)## [1] 128.9299Plot the model
bs2506_SLW_lm_plottemp <- ggplot(bsdf_2506_SL_W, aes(x = SL_CM, y = W_G)) +
geom_point(aes(SL_CM, W_G),
color = "#242124",
alpha = 0.8) +
labs(x = "Fork Length in cm",
y = "Weight in grams",
title = "Blackfin Snapper (n = 33)") +
theme_classic() +
geom_smooth(method = lm,
aes(SL_CM, W_G),
linewidth = 0.5,
se = FALSE,
color = safe_colorblind_palette[3])
# add label with equation
bs2506_SLW_lm_plot <- bs2506_SLW_lm_plottemp +
geom_text(x = 16, y = 1400,
label = expression(W~'(g)' == 0.020 * TL~'(cm)' + 21.479~"," ~ R^{2} == 0.935 ~ "(Linear)"),
colour = safe_colorblind_palette[3],
hjust = 0)
bs2506_SLW_lm_plot
NLS Models
NLS 1:
bs2506_SLW_nls1 <- nls_multstart(W_G~A*SL_CM^B,
data = bsdf_2506_SL_W,
iter = 500,
start_lower = c(A=0, B=2),
start_upper = c(A=0.5, B=3.5))
summary(bs2506_SLW_nls1)##
## Formula: W_G ~ A * SL_CM^B
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## A 0.016531 0.006721 2.46 0.0197 *
## B 3.126597 0.115643 27.04 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 66.22 on 31 degrees of freedom
##
## Number of iterations to convergence: 9
## Achieved convergence tolerance: 1.49e-08Plot a histogram of the residuals to check normality of residuals
qqnorm(resid(bs2506_SLW_nls1))
qqline(resid(bs2506_SLW_nls1))
Heavy tailed but okay
Calculate R-squared
rss_fitted <- sum(residuals(bs2506_SLW_nls1)^2)
rss_constant <- sum((bsdf_2506_SL_W$W_G - mean(bsdf_2506_SL_W$W_G))^2)
approx_r_squared <- 1 - rss_fitted / rss_constant
approx_r_squared## [1] 0.9726183Calculate AIC
AIC(bs2506_SLW_nls1)## [1] 374.3251Calculate the predicted W_kg values based on this NLS model
bs2506_SLW_nls1_x <- c(16:39) #rounded actual min and max SL_CM values
bs2506_SLW_nls1_predic <- data.frame(bs2506_SLW_nls1_x)
bs2506_SLW_nls1_f <- function(bs2506_SLW_nls1_x) {
#formula using calculated parameters
return(0.016531*bs2506_SLW_nls1_x^3.126597)
}
bs2506_SLW_nls1_predic$bs2506_SLW_nls1_y <- sapply(bs2506_SLW_nls1_predic$bs2506_SLW_nls1_x, bs2506_SLW_nls1_f)
View(bs2506_SLW_nls1_predic)Plot the model
bs2506_SLW_nls1_plot_title <- expression(italic("L. buccanella")~"SL-W Relationship (n = 33)")
bs2506_SLW_nls1_plot_temp <- ggplot(bsdf_2506_SL_W, aes(x = SL_CM, y = W_G)) +
geom_point(aes(SL_CM, W_G),
color = "#242124",
alpha = 0.8) +
labs(x = expression("Standard Length in cm"),
y = expression("Weight in grams"),
title = bs2506_SLW_nls1_plot_title) +
geom_smooth(aes(x = bs2506_SLW_nls1_x, y = bs2506_SLW_nls1_y),
data = bs2506_SLW_nls1_predic,
method = "loess",
se = F,
color = "#73dfeb",
linewidth = 0.75) +
theme_classic()
bs2506_SLW_nls1_plot <- bs2506_SLW_nls1_plot_temp +
geom_text(x = 16, y = 1500,
label = expression(W~'(g)' == 1.6531%*%10^{-02} * TL~'(cm)'^{3.126}~"," ~ R^{2} == 0.973 ~ '(NLS)'),
colour = "#73dfeb",
hjust = 0)
bs2506_SLW_nls1_plot
Maturity
Load data file
bsdf_2506_maturity <- read.csv("bs2506df.csv")
bsdf_2506_maturity$catch_month <- factor(bsdf_2506_maturity$catch_month)
bsdf_2506_maturity$year <- factor(bsdf_2506_maturity$year)
bsdf_2506_maturity$species <- factor(bsdf_2506_maturity$species)
bsdf_2506_maturity$common_name <- factor(bsdf_2506_maturity$common_name)
bsdf_2506_maturity$frozen <- factor(bsdf_2506_maturity$frozen)
bsdf_2506_maturity$catch_method <- factor(bsdf_2506_maturity$catch_method)
bsdf_2506_maturity$bait <- factor(bsdf_2506_maturity$bait)
bsdf_2506_maturity$area <- factor(bsdf_2506_maturity$area)
bsdf_2506_maturity$gutted <- factor(bsdf_2506_maturity$gutted)
bsdf_2506_maturity$sex <- factor(bsdf_2506_maturity$sex)
bsdf_2506_maturity$maturity <- factor(bsdf_2506_maturity$maturity)
#drop the columns with neither sex nor maturity measurements
bsdf_2506_maturity <- bsdf_2506_maturity %>%
filter(maturity != "" | sex != "")
View(bsdf_2506_maturity)
summary(bsdf_2506_maturity)## X catch_date catch_month landing_time
## Min. : 1.00 Length:55 : 2 Length:55
## 1st Qu.:14.50 Class :character Apr:11 Class :character
## Median :28.00 Mode :character Dec: 1 Mode :character
## Mean :36.93 Jul: 6
## 3rd Qu.:66.00 Jun: 8
## Max. :85.00 May:19
## Sep: 8
## dissection_date dissection_time year month
## Length:55 Length:55 2023:22 Length:55
## Class :character Class :character 2024:33 Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
## observer species common_name
## Length:55 Lutjanus buccanella:55 Blackfin Snapper:55
## Class :character
## Mode :character
##
##
##
##
## fish_code frozen catch_method bait
## Length:55 N:33 :21 :21
## Class :character Y:22 Lobster trap:11 Cowskin :11
## Mode :character Redfish trap:23 Japanese bait :22
## Japanese bait; Squid: 1
##
##
##
## area SL_CM FL_CM TL_CM HH_CM
## A3/A4 :17 Min. :16.40 Min. :18.70 Min. :20.50 Min. :3.600
## :13 1st Qu.:24.60 1st Qu.:24.20 1st Qu.:25.60 1st Qu.:6.200
## D3 :10 Median :27.95 Median :27.90 Median :29.80 Median :7.950
## C2 : 8 Mean :28.33 Mean :29.28 Mean :31.48 Mean :7.359
## D5 : 3 3rd Qu.:34.08 3rd Qu.:32.85 3rd Qu.:35.50 3rd Qu.:8.850
## D4 : 2 Max. :38.30 Max. :42.90 Max. :45.90 Max. :9.900
## (Other): 2 NA's :21 NA's :2 NA's :33
## ED_CM W_G finclip otoliths gutted
## Min. :1.200 Min. : 108.5 Length:55 Length:55 N:54
## 1st Qu.:2.200 1st Qu.: 258.1 Class :character Class :character Y: 1
## Median :2.400 Median : 364.6 Mode :character Mode :character
## Mean :2.245 Mean : 514.6
## 3rd Qu.:2.500 3rd Qu.: 638.8
## Max. :2.900 Max. :1493.9
## NA's :33 NA's :2
## gonads_present sex maturity G1L_CM G2L_CM
## Length:55 : 2 :22 Min. :1.600 Min. :1.200
## Class :character F:19 A:15 1st Qu.:2.625 1st Qu.:2.700
## Mode :character M:34 B:10 Median :3.500 Median :3.600
## C: 8 Mean :3.915 Mean :3.908
## 3rd Qu.:4.500 3rd Qu.:4.600
## Max. :9.300 Max. :8.800
## NA's :1 NA's :6
## GW_G stomach._sample
## Length:55 Length:55
## Class :character Class :character
## Mode :character Mode :character
##
##
##
## calculate sex ratio
bs_totalfemale <- bsdf_2506_maturity %>%
filter(sex == "F") %>%
count()
bs_totalfemale## n
## 1 19bs_totalmale <- bsdf_2506_maturity %>%
filter(sex == "M") %>%
count()
bs_totalmale## n
## 1 34bs_sexratio <- bs_totalfemale/bs_totalmale
bs_sexratio #female number## n
## 1 0.5588235Descriptive statistics
bs_sex_summary <- bsdf_2506_maturity %>%
group_by(sex) %>%
summarize(
count = n(),
mean = mean(FL_CM, na.rm = TRUE),
sd = sd(FL_CM, na.rm = TRUE),
median = median(FL_CM, na.rm = TRUE),
min = min(FL_CM, na.rm = TRUE),
max = max(FL_CM, na.rm = TRUE),
IQR = IQR(FL_CM, na.rm = TRUE)
)
print(bs_sex_summary, digits = 4)## # A tibble: 3 × 8
## sex count mean sd median min max IQR
## <fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 "" 2 24.6 8.41 24.6 18.7 30.6 5.95
## 2 "F" 19 25.4 3.71 25.9 19.4 34.3 5.1
## 3 "M" 34 31.7 6.69 30.7 20.8 42.9 11.4bsdf_2506_maturity <- bsdf_2506_maturity %>%
mutate(maturity = as.character(maturity)) %>%
mutate(maturity_state = case_when(maturity == "A" ~ "immature",
maturity == "B" ~ "mature",
maturity == "C" ~ "mature",
maturity == "D" ~ "mature")) %>%
mutate(maturity = as.factor(maturity),
maturity_state = as.factor(maturity_state))
bs_maturity_summary <- bsdf_2506_maturity %>%
group_by(maturity_state) %>%
summarize(
count = n(),
mean = mean(FL_CM, na.rm = TRUE),
sd = sd(FL_CM, na.rm = TRUE),
median = median(FL_CM, na.rm = TRUE),
min = min(FL_CM, na.rm = TRUE),
max = max(FL_CM, na.rm = TRUE),
IQR = IQR(FL_CM, na.rm = TRUE)
)
bs_maturity_summary## # A tibble: 3 × 8
## maturity_state count mean sd median min max IQR
## <fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 immature 15 28.3 4.29 28.6 18.7 34.3 3.45
## 2 mature 18 35.3 6.64 38.7 23.9 42.9 12.0
## 3 <NA> 22 25.1 3.53 24 20.5 33.3 4.62