L. vivanus body size models
Naomi Bakken
2024-08-15
This R markdown file includes all the R code used to analyze Yelloweye Snapper (Lutjanus vivanus) 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
ysdf_2506_LLrelationships <- read.csv("ys2506df.csv")
ysdf_2506_LLrelationships$catch_month <- factor(ysdf_2506_LLrelationships$catch_month)
ysdf_2506_LLrelationships$year <- factor(ysdf_2506_LLrelationships$year)
ysdf_2506_LLrelationships$species <- factor(ysdf_2506_LLrelationships$species)
ysdf_2506_LLrelationships$common_name <- factor(ysdf_2506_LLrelationships$common_name)
ysdf_2506_LLrelationships$frozen <- factor(ysdf_2506_LLrelationships$frozen)
ysdf_2506_LLrelationships$catch_method <- factor(ysdf_2506_LLrelationships$catch_method)
ysdf_2506_LLrelationships$bait <- factor(ysdf_2506_LLrelationships$bait)
ysdf_2506_LLrelationships$area <- factor(ysdf_2506_LLrelationships$area)
ysdf_2506_LLrelationships$gutted <- factor(ysdf_2506_LLrelationships$gutted)
ysdf_2506_LLrelationships$sex <- factor(ysdf_2506_LLrelationships$sex)
ysdf_2506_LLrelationships$maturity <- factor(ysdf_2506_LLrelationships$maturity)
summary(ysdf_2506_LLrelationships)## X catch_date catch_month landing_time
## Min. : 1.00 Length:86 Sep :32 Length:86
## 1st Qu.:22.25 Class :character Jun :19 Class :character
## Median :43.50 Mode :character May :14 Mode :character
## Mean :43.50 Jul : 9
## 3rd Qu.:64.75 Mar : 5
## Max. :86.00 : 3
## (Other): 4
## dissection_date dissection_time year month
## Length:86 Length:86 2023:57 Length:86
## Class :character Class :character 2024:29 Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
## observer species common_name
## Length:86 Lutjanus vivanus:86 Yelloweye Snapper:86
## Class :character
## Mode :character
##
##
##
##
## fish_code frozen catch_method bait
## Length:86 N:26 :48 :54
## Class :character Y:60 Dropline : 8 Cowskin : 1
## Mode :character Lobster trap: 1 Japanese bait :23
## Redfish trap:29 Japanese bait; Squid: 8
##
##
##
## area SL_CM FL_CM TL_CM HH_CM
## :25 Min. :16.00 Min. :17.20 Min. :18.30 Min. :3.600
## A3/A4: 5 1st Qu.:24.68 1st Qu.:21.73 1st Qu.:23.20 1st Qu.:5.250
## B4 :13 Median :29.20 Median :26.15 Median :28.60 Median :6.600
## C2 :24 Mean :28.70 Mean :27.57 Mean :29.77 Mean :6.317
## D3 : 6 3rd Qu.:32.10 3rd Qu.:33.08 3rd Qu.:36.08 3rd Qu.:7.450
## D5 :13 Max. :40.40 Max. :45.90 Max. :49.20 Max. :9.400
## NA's :48 NA's :63
## ED_CM W_G finclip otoliths gutted
## Min. :1.20 Min. : 83.9 Length:86 Length:86 N:86
## 1st Qu.:1.65 1st Qu.: 162.7 Class :character Class :character
## Median :1.90 Median : 281.9 Mode :character Mode :character
## Mean :1.87 Mean : 429.0
## 3rd Qu.:2.10 3rd Qu.: 583.3
## Max. :3.00 Max. :1714.7
## NA's :63
## gonads_present sex maturity G1L_CM G2L_CM
## Length:86 :43 :57 Min. : 1.000 Min. : 0.000
## Class :character F:22 A:13 1st Qu.: 2.250 1st Qu.: 2.200
## Mode :character M:21 B:12 Median : 3.800 Median : 3.600
## C: 2 Mean : 4.441 Mean : 4.261
## D: 2 3rd Qu.: 5.950 3rd Qu.: 5.500
## Max. :11.500 Max. :11.600
## NA's :27 NA's :29
## GW_G stomach._sample
## Min. :0.0000 Length:86
## 1st Qu.:0.0000 Class :character
## Median :0.4000 Mode :character
## Mean :0.9695
## 3rd Qu.:1.2500
## Max. :5.6000
## NA's :27View(ysdf_2506_LLrelationships)
sumtable(ysdf_2506_LLrelationships, digits = 4)| Variable | N | Mean | Std. Dev. | Min | Pctl. 25 | Pctl. 75 | Max |
|---|---|---|---|---|---|---|---|
| X | 86 | 43.5 | 24.97 | 1 | 22.25 | 64.75 | 86 |
| catch_month | 86 | ||||||
| … | 3 | 3.49% | |||||
| … Apr | 1 | 1.16% | |||||
| … Aug | 2 | 2.33% | |||||
| … Dec | 1 | 1.16% | |||||
| … Jul | 9 | 10.47% | |||||
| … Jun | 19 | 22.09% | |||||
| … Mar | 5 | 5.81% | |||||
| … May | 14 | 16.28% | |||||
| … Sep | 32 | 37.21% | |||||
| year | 86 | ||||||
| … 2023 | 57 | 66.28% | |||||
| … 2024 | 29 | 33.72% | |||||
| species | 86 | ||||||
| … Lutjanus vivanus | 86 | 100% | |||||
| common_name | 86 | ||||||
| … Yelloweye Snapper | 86 | 100% | |||||
| frozen | 86 | ||||||
| … N | 26 | 30.23% | |||||
| … Y | 60 | 69.77% | |||||
| catch_method | 86 | ||||||
| … | 48 | 55.81% | |||||
| … Dropline | 8 | 9.3% | |||||
| … Lobster trap | 1 | 1.16% | |||||
| … Redfish trap | 29 | 33.72% | |||||
| bait | 86 | ||||||
| … | 54 | 62.79% | |||||
| … Cowskin | 1 | 1.16% | |||||
| … Japanese bait | 23 | 26.74% | |||||
| … Japanese bait; Squid | 8 | 9.3% | |||||
| area | 86 | ||||||
| … | 25 | 29.07% | |||||
| … A3/A4 | 5 | 5.81% | |||||
| … B4 | 13 | 15.12% | |||||
| … C2 | 24 | 27.91% | |||||
| … D3 | 6 | 6.98% | |||||
| … D5 | 13 | 15.12% | |||||
| SL_CM | 38 | 28.7 | 6.438 | 16 | 24.68 | 32.1 | 40.4 |
| FL_CM | 86 | 27.57 | 7.362 | 17.2 | 21.72 | 33.08 | 45.9 |
| TL_CM | 86 | 29.77 | 7.916 | 18.3 | 23.2 | 36.08 | 49.2 |
| HH_CM | 23 | 6.317 | 1.622 | 3.6 | 5.25 | 7.45 | 9.4 |
| ED_CM | 23 | 1.87 | 0.4311 | 1.2 | 1.65 | 2.1 | 3 |
| W_G | 86 | 429 | 390.9 | 83.9 | 162.7 | 583.3 | 1715 |
| finclip | 86 | ||||||
| … | 47 | 54.65% | |||||
| … N | 5 | 5.81% | |||||
| … Y | 34 | 39.53% | |||||
| otoliths | 86 | ||||||
| … y | 1 | 1.16% | |||||
| … Y | 84 | 97.67% | |||||
| … Y | 1 | 1.16% | |||||
| gutted | 86 | ||||||
| … N | 86 | 100% | |||||
| gonads_present | 86 | ||||||
| … N | 28 | 32.56% | |||||
| … y | 1 | 1.16% | |||||
| … Y | 57 | 66.28% | |||||
| sex | 86 | ||||||
| … | 43 | 50% | |||||
| … F | 22 | 25.58% | |||||
| … M | 21 | 24.42% | |||||
| maturity | 86 | ||||||
| … | 57 | 66.28% | |||||
| … A | 13 | 15.12% | |||||
| … B | 12 | 13.95% | |||||
| … C | 2 | 2.33% | |||||
| … D | 2 | 2.33% | |||||
| G1L_CM | 59 | 4.441 | 2.899 | 1 | 2.25 | 5.95 | 11.5 |
| G2L_CM | 57 | 4.261 | 2.963 | 0 | 2.2 | 5.5 | 11.6 |
| GW_G | 59 | 0.9695 | 1.377 | 0 | 0 | 1.25 | 5.6 |
| stomach._sample | 86 | ||||||
| … N | 62 | 72.09% | |||||
| … Y | 24 | 27.91% |
Fork length (FL_CM) - Total length (TL_CM) relationship
ys2506_FLTL_plotbase <- ggplot(ysdf_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 vivanus FL-TL relationship (n = 86)",
fill = "Sex") +
scale_fill_viridis(option = "C",
discrete = TRUE,
begin = 0.2,
end = 0.8,
labels = c("N/A", "Female", "Male"))
ys2506_FLTL_plotbase
Linear model
ys2506_FLTL_lm <- lm(FL_CM ~ TL_CM, data = ysdf_2506_LLrelationships)
summary(ys2506_FLTL_lm)##
## Call:
## lm(formula = FL_CM ~ TL_CM, data = ysdf_2506_LLrelationships)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.1071 -0.2900 -0.0765 0.1598 8.5843
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.15871 0.44125 0.36 0.72
## TL_CM 0.92076 0.01433 64.25 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.046 on 84 degrees of freedom
## Multiple R-squared: 0.9801, Adjusted R-squared: 0.9798
## F-statistic: 4127 on 1 and 84 DF, p-value: < 2.2e-16length(ys2506_FLTL_lm$fitted.values)## [1] 86Scatterplot the residuals against the fitted values to check homogeneity of variance
plot(resid(ys2506_FLTL_lm) ~ fitted(ys2506_FLTL_lm))
Fine other than an outlier
Plot a histogram to check normality of residuals
qqnorm(ys2506_FLTL_lm$residuals)
qqline(ys2506_FLTL_lm$residuals)
fine
Fork length (FL_CM) - Standard length (SL_CM) relationship
ys2506_FLSL_plot <- ggplot(ysdf_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 vivanus FL-SL relationship (n = 38)",
fill = "Sex") +
scale_fill_viridis(option = "C",
discrete = TRUE,
begin = 0.2,
end = 0.8,
labels = c("N/A", "Female", "Male"))
ys2506_FLSL_plot
Linear model
ys2506_FLSL_lm <- lm(FL_CM ~ SL_CM, data = ysdf_2506_LLrelationships)
print(summary(ys2506_FLSL_lm), digits = 5)##
## Call:
## lm(formula = FL_CM ~ SL_CM, data = ysdf_2506_LLrelationships)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.289805 -0.319391 0.030364 0.410727 0.956810
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.231617 0.446676 0.5185 0.6073
## SL_CM 1.138100 0.015197 74.8892 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.59514 on 36 degrees of freedom
## (48 observations deleted due to missingness)
## Multiple R-squared: 0.99362, Adjusted R-squared: 0.99344
## F-statistic: 5608.4 on 1 and 36 DF, p-value: < 2.22e-16length(ys2506_FLSL_lm$fitted.values)## [1] 38Scatterplot the residuals against the fitted values to check homogeneity of variance
plot(resid(ys2506_FLSL_lm) ~ fitted(ys2506_FLSL_lm))
Plot a histogram to check normality of residuals
qqnorm(resid(ys2506_FLSL_lm))
qqline(resid(ys2506_FLSL_lm))
Total length (TL_CM) - Standard length (SL_CM) relationship
ys2506_TLSL_plot <- ggplot(ysdf_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 vivanus TL-SL relationship (n = 38)",
fill = "Sex") +
scale_fill_viridis(option = "C",
discrete = TRUE,
begin = 0.2,
end = 0.8,
labels = c("N/A", "Female", "Male"))
ys2506_FLSL_plot
Linear model
ys2506_TLSL_lm <- lm(TL_CM ~ SL_CM, data = ysdf_2506_LLrelationships)
summary(ys2506_TLSL_lm)##
## Call:
## lm(formula = TL_CM ~ SL_CM, data = ysdf_2506_LLrelationships)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.7598 -0.1452 0.1529 0.8277 1.4922
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.22840 1.28058 0.959 0.344
## SL_CM 1.18993 0.04357 27.312 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.706 on 36 degrees of freedom
## (48 observations deleted due to missingness)
## Multiple R-squared: 0.954, Adjusted R-squared: 0.9527
## F-statistic: 745.9 on 1 and 36 DF, p-value: < 2.2e-16length(ys2506_TLSL_lm$fitted.values)## [1] 38Scatterplot the residuals against the fitted values to check homogeneity of variance
plot(resid(ys2506_TLSL_lm) ~ fitted(ys2506_TLSL_lm))
Plot a histogram to check normality of residuals
qqnorm(resid(ys2506_TLSL_lm))
qqline(resid(ys2506_TLSL_lm))
Fork length (FL_CM) - Weight (W_G) relationship
Load data file
ysdf_2506_FL_W <- read.csv("ys2506df.csv")
ysdf_2506_FL_W$catch_month <- factor(ysdf_2506_FL_W$catch_month)
ysdf_2506_FL_W$year <- factor(ysdf_2506_FL_W$year)
ysdf_2506_FL_W$species <- factor(ysdf_2506_FL_W$species)
ysdf_2506_FL_W$common_name <- factor(ysdf_2506_FL_W$common_name)
ysdf_2506_FL_W$frozen <- factor(ysdf_2506_FL_W$frozen)
ysdf_2506_FL_W$catch_method <- factor(ysdf_2506_FL_W$catch_method)
ysdf_2506_FL_W$bait <- factor(ysdf_2506_FL_W$bait)
ysdf_2506_FL_W$area <- factor(ysdf_2506_FL_W$area)
ysdf_2506_FL_W$gutted <- factor(ysdf_2506_FL_W$gutted)
ysdf_2506_FL_W$sex <- factor(ysdf_2506_FL_W$sex)
ysdf_2506_FL_W$maturity <- factor(ysdf_2506_FL_W$maturity)
# Remove gutted samples
ysdf_2506_FL_W <- subset(filter(ysdf_2506_FL_W, gutted == "N"))
# remove NA rows
ysdf_2506_FL_W <- ysdf_2506_FL_W%>%
drop_na(FL_CM)
View(ysdf_2506_FL_W)
summary(ysdf_2506_FL_W)## X catch_date catch_month landing_time
## Min. : 1.00 Length:86 Sep :32 Length:86
## 1st Qu.:22.25 Class :character Jun :19 Class :character
## Median :43.50 Mode :character May :14 Mode :character
## Mean :43.50 Jul : 9
## 3rd Qu.:64.75 Mar : 5
## Max. :86.00 : 3
## (Other): 4
## dissection_date dissection_time year month
## Length:86 Length:86 2023:57 Length:86
## Class :character Class :character 2024:29 Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
## observer species common_name
## Length:86 Lutjanus vivanus:86 Yelloweye Snapper:86
## Class :character
## Mode :character
##
##
##
##
## fish_code frozen catch_method bait
## Length:86 N:26 :48 :54
## Class :character Y:60 Dropline : 8 Cowskin : 1
## Mode :character Lobster trap: 1 Japanese bait :23
## Redfish trap:29 Japanese bait; Squid: 8
##
##
##
## area SL_CM FL_CM TL_CM HH_CM
## :25 Min. :16.00 Min. :17.20 Min. :18.30 Min. :3.600
## A3/A4: 5 1st Qu.:24.68 1st Qu.:21.73 1st Qu.:23.20 1st Qu.:5.250
## B4 :13 Median :29.20 Median :26.15 Median :28.60 Median :6.600
## C2 :24 Mean :28.70 Mean :27.57 Mean :29.77 Mean :6.317
## D3 : 6 3rd Qu.:32.10 3rd Qu.:33.08 3rd Qu.:36.08 3rd Qu.:7.450
## D5 :13 Max. :40.40 Max. :45.90 Max. :49.20 Max. :9.400
## NA's :48 NA's :63
## ED_CM W_G finclip otoliths gutted
## Min. :1.20 Min. : 83.9 Length:86 Length:86 N:86
## 1st Qu.:1.65 1st Qu.: 162.7 Class :character Class :character
## Median :1.90 Median : 281.9 Mode :character Mode :character
## Mean :1.87 Mean : 429.0
## 3rd Qu.:2.10 3rd Qu.: 583.3
## Max. :3.00 Max. :1714.7
## NA's :63
## gonads_present sex maturity G1L_CM G2L_CM
## Length:86 :43 :57 Min. : 1.000 Min. : 0.000
## Class :character F:22 A:13 1st Qu.: 2.250 1st Qu.: 2.200
## Mode :character M:21 B:12 Median : 3.800 Median : 3.600
## C: 2 Mean : 4.441 Mean : 4.261
## D: 2 3rd Qu.: 5.950 3rd Qu.: 5.500
## Max. :11.500 Max. :11.600
## NA's :27 NA's :29
## GW_G stomach._sample
## Min. :0.0000 Length:86
## 1st Qu.:0.0000 Class :character
## Median :0.4000 Mode :character
## Mean :0.9695
## 3rd Qu.:1.2500
## Max. :5.6000
## NA's :27First look at the data with a scatterplot
fill = sex
ggplot(ysdf_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 vivanus, n = 86") +
scale_fill_viridis(option = "D", discrete = TRUE)
Check for outliers in data: calculate z-scores
yszscore_FL_CM <- (ysdf_2506_FL_W$FL_CM - mean(ysdf_2506_FL_W$FL_CM)) / sd(ysdf_2506_FL_W$FL_CM)
ysdf_2506_FL_W$zscore_FL_CM <- yszscore_FL_CM
yszscore_Wgrams <- (ysdf_2506_FL_W$W_G - mean(ysdf_2506_FL_W$W_G)) / sd(ysdf_2506_FL_W$W_G)
ysdf_2506_FL_W$zscore_Wgrams <- yszscore_Wgrams
View(ysdf_2506_FL_W)z-scores above 3 are generally considered to be outliers -> remove from the dataset
Filter the data
ysdf_2506_FL_W <- ysdf_2506_FL_W[ysdf_2506_FL_W$zscore_Wgrams < 3, ]
View(ysdf_2506_FL_W)
summary(ysdf_2506_FL_W)## X catch_date catch_month landing_time
## Min. : 1.00 Length:84 Sep :32 Length:84
## 1st Qu.:23.75 Class :character Jun :18 Class :character
## Median :44.50 Mode :character May :13 Mode :character
## Mean :44.13 Jul : 9
## 3rd Qu.:65.25 Mar : 5
## Max. :86.00 : 3
## (Other): 4
## dissection_date dissection_time year month
## Length:84 Length:84 2023:57 Length:84
## Class :character Class :character 2024:27 Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
## observer species common_name
## Length:84 Lutjanus vivanus:84 Yelloweye Snapper:84
## Class :character
## Mode :character
##
##
##
##
## fish_code frozen catch_method bait
## Length:84 N:24 :48 :54
## Class :character Y:60 Dropline : 8 Cowskin : 1
## Mode :character Lobster trap: 1 Japanese bait :21
## Redfish trap:27 Japanese bait; Squid: 8
##
##
##
## area SL_CM FL_CM TL_CM HH_CM
## :25 Min. :16.00 Min. :17.20 Min. :18.30 Min. :3.600
## A3/A4: 5 1st Qu.:24.32 1st Qu.:21.65 1st Qu.:23.18 1st Qu.:4.900
## B4 :13 Median :28.70 Median :25.55 Median :27.95 Median :6.500
## C2 :24 Mean :28.24 Mean :27.23 Mean :29.51 Mean :6.124
## D3 : 6 3rd Qu.:31.65 3rd Qu.:32.62 3rd Qu.:35.55 3rd Qu.:6.700
## D5 :11 Max. :40.40 Max. :45.90 Max. :49.20 Max. :9.400
## NA's :48 NA's :63
## ED_CM W_G finclip otoliths gutted
## Min. :1.20 Min. : 83.9 Length:84 Length:84 N:84
## 1st Qu.:1.60 1st Qu.: 161.3 Class :character Class :character
## Median :1.90 Median : 272.4 Mode :character Mode :character
## Mean :1.81 Mean : 399.6
## 3rd Qu.:2.10 3rd Qu.: 562.1
## Max. :2.50 Max. :1553.3
## NA's :63
## gonads_present sex maturity G1L_CM G2L_CM
## Length:84 :43 :57 Min. : 1.000 Min. : 0.000
## Class :character F:21 A:12 1st Qu.: 2.200 1st Qu.: 2.050
## Mode :character M:20 B:11 Median : 3.700 Median : 3.200
## C: 2 Mean : 4.195 Mean : 4.002
## D: 2 3rd Qu.: 5.700 3rd Qu.: 5.300
## Max. :10.900 Max. :11.100
## NA's :27 NA's :29
## GW_G stomach._sample zscore_FL_CM zscore_Wgrams
## Min. :0.0000 Length:84 Min. :-1.40788 Min. :-0.88288
## 1st Qu.:0.0000 Class :character 1st Qu.:-0.80344 1st Qu.:-0.68482
## Median :0.3000 Mode :character Median :-0.27371 Median :-0.40054
## Mean :0.8807 Mean :-0.04523 Mean :-0.07534
## 3rd Qu.:1.1000 3rd Qu.: 0.68728 3rd Qu.: 0.34036
## Max. :5.6000 Max. : 2.49040 Max. : 2.87605
## NA's :27LM Linear model FL-W
ys2506_FLW_lm <- lm(FL_CM ~ W_G, data = ysdf_2506_FL_W)
summary(ys2506_FLW_lm)##
## Call:
## lm(formula = FL_CM ~ W_G, data = ysdf_2506_FL_W)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.0801 -1.4795 0.2457 1.5304 3.1382
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.932e+01 3.365e-01 57.41 <2e-16 ***
## W_G 1.980e-02 6.394e-04 30.97 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.007 on 82 degrees of freedom
## Multiple R-squared: 0.9212, Adjusted R-squared: 0.9203
## F-statistic: 959.1 on 1 and 82 DF, p-value: < 2.2e-16Scatterplot the residuals against the fitted values to check homogeneity of variance
plot(resid(ys2506_FLW_lm) ~ fitted(ys2506_FLW_lm))
The scatter plot clearly looks like it indicates a nonlinear
relationship
Plot a histogram to check normality of residuals
qqnorm(resid(ys2506_FLW_lm))
qqline(resid(ys2506_FLW_lm))
Plot the model
ys2506_FLW_lm_plottemp <- ggplot(ysdf_2506_FL_W, aes(x = FL_CM, y = W_G)) +
geom_point(aes(FL_CM, W_G),
color = "black",
alpha = 0.6) +
labs(x = "Fork Length in cm",
y = "Weight in grams",
title = "Lutjanus vivanus, n = 84") +
theme_classic() +
geom_smooth(method = lm,
aes(FL_CM, W_G),
linewidth = 0.5,
se = FALSE,
color = safe_colorblind_palette[1])
# add label with equation
ys2506_FLW_lm_plot <- ys2506_FLW_lm_plottemp +
geom_text(x = 18, y = 1500,
label = expression(W~'(g)' == 0.031 * FL~'(cm)' + 16.444~"," ~ r^{2} == 0.954 ~ "(Linear)"),
colour = safe_colorblind_palette[1],
hjust = 0)
ys2506_FLW_lm_plot
NLS Models
NLS 1
ys2506_FLW_nls <- nls_multstart(W_G~A*FL_CM^B,
data = ysdf_2506_FL_W,
iter = 500,
start_lower = c(A=0, B=2),
start_upper = c(A=0.5, B=3.5))
summary(ys2506_FLW_nls)##
## Formula: W_G ~ A * FL_CM^B
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## A 0.013323 0.002134 6.244 1.79e-08 ***
## B 3.057844 0.044210 69.167 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 40.62 on 82 degrees of freedom
##
## Number of iterations to convergence: 17
## Achieved convergence tolerance: 1.49e-08Plot a histogram of the residuals to check normality of residuals
#check residuals
qqnorm(resid(ys2506_FLW_nls))
qqline(resid(ys2506_FLW_nls))
heavy tails but okay
Calculate R-squared
rss_fitted <- sum(residuals(ys2506_FLW_nls)^2)
rss_constant <- sum((ysdf_2506_FL_W$W_G - mean(ysdf_2506_FL_W$W_G))^2)
approx_r_squared <- 1 - rss_fitted / rss_constant
approx_r_squared## [1] 0.9862662Calculate AIC
AIC(ys2506_FLW_nls)## [1] 864.6691Calculate the predicted W_kg values based on this NLS model
ys2506_FLW_nls_x <- c(17:46) #rounded actual min and max FL_MM values
ys2506_FLW_nls_predic <- data.frame(ys2506_FLW_nls_x)
ys2506_FLW_nls_f <- function(ys2506_FLW_nls_x) {
#formula using calculated parameters
return(0.013323*ys2506_FLW_nls_x^3.057844)
}
ys2506_FLW_nls_predic$ys2506_FLW_nls_y <- sapply(ys2506_FLW_nls_predic$ys2506_FLW_nls_x, ys2506_FLW_nls_f)
View(ys2506_FLW_nls_predic)Plot the model on the same graph as the other model
ys2506_FLW_nls_plottemp <- ys2506_FLW_lm_plot +
geom_smooth(aes(x = ys2506_FLW_nls_x, y = ys2506_FLW_nls_y),
data = ys2506_FLW_nls_predic,
method = "loess",
se = F,
color = safe_colorblind_palette[2],
linewidth = 0.75,
alpha = 0.75)
# add label with equation
ys2506_FLW_nls_plot <- ys2506_FLW_nls_plottemp +
geom_text(x = 18, y = 1400,
label = expression(W~'(g)' == 1.332%*%10^{-02} * FL~'(cm)'^{3.058}~"," ~ R^{2} == 0.986 ~ '(NLS)'),
colour = safe_colorblind_palette[2],
hjust = 0)
ys2506_FLW_nls_plot
Brms1: with nls parameters and a wide normal distribution around them as priors
set priors and run model
#set priors
ys2506_FLW_prior1 <- prior(normal(0.013323, 0.005) ,nlpar = "A")+prior(normal(3.057844, 1),nlpar = "B")
ys2506_FLW_brm1 <- brm(bf(W_G~A*FL_CM^B, A + B ~ 1, nl = TRUE),data = ysdf_2506_FL_W,
prior = ys2506_FLW_prior1,
warmup = 1000,
iter = 5000,
chains = 4,
cores = 4,
control = list(adapt_delta = 0.8))
print(summary(ys2506_FLW_brm1), digits = 6)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: W_G ~ A * FL_CM^B
## A ~ 1
## B ~ 1
## Data: ysdf_2506_FL_W (Number of observations: 84)
## Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 1;
## total post-warmup draws = 16000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## A_Intercept 0.013743 0.002019 0.010132 0.018013 1.000656 4764 4774
## B_Intercept 3.052201 0.040566 2.974733 3.133695 1.000494 4757 4780
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 41.200042 3.331157 35.306745 48.355494 1.000616 5461 5942
##
## 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).Plot a histogram of the residuals to check normality of residuals
qqnorm(resid(ys2506_FLW_brm1))
qqline(resid(ys2506_FLW_brm1))
normal enough
Calculate R-squared
rss_fitted <- sum(residuals(ys2506_FLW_brm1)^2)
rss_constant <- sum((ysdf_2506_FL_W$W_G - mean(ysdf_2506_FL_W$W_G))^2)
approx_r_squared <- 1 - rss_fitted / rss_constant
approx_r_squared ## [1] 0.8289979smaller than the nls r-squared
Calculate the predicted W_G values based on this NLS model
ys2506_FLW_brm1_x <- c(17:46) #rounded actual min and max TL values
ys2506_FLW_brm1_predic <- data.frame(ys2506_FLW_brm1_x)
ys2506_FLW_brm1_f <- function(ys2506_FLW_brm1_x) {
#formula using calculated parameters
return(0.013688*ys2506_FLW_brm1_x^3.053169)
}
ys2506_FLW_brm1_predic$ys2506_FLW_brm1_y <- sapply(ys2506_FLW_brm1_predic$ys2506_FLW_brm1_x, ys2506_FLW_brm1_f)
View(ys2506_FLW_brm1_predic)Plot
ys2506_FLW_brm1_plottemp <- ys2506_FLW_nls_plot +
geom_smooth(aes(x = ys2506_FLW_brm1_x, y = ys2506_FLW_brm1_y),
data = ys2506_FLW_brm1_predic,
method = "loess",
se = F,
color = safe_colorblind_palette[4],
linewidth = 0.75,
alpha = 0.75)
# add label with equation
ys2506_FLW_brm1_plot <- ys2506_FLW_brm1_plottemp +
geom_text(x = 18, y = 1300,
label = expression(W~'(g)' == 1.369%*%10^{-02} * FL~'(cm)'^{3.053}~"," ~ R^{2} == 0.829 ~ '(BRMS)'),
colour = safe_colorblind_palette[4],
hjust = 0)
ys2506_FLW_brm1_plot
Brms2: with previous brms parameters as priors
also tightening the normal distribution around the starting priors using the difference between the upper and lower CI from the previous models
#A UCI-LCI
0.017918-0.010158 ## [1] 0.00776#B UCI-LCI
3.132711-2.976053## [1] 0.156658set priors and run model
#set priors
ys2506_FLW_prior2 <- prior(normal(0.013688, 0.00776) ,nlpar = "A")+prior(normal(3.053169, 0.156658),nlpar = "B")
ys2506_FLW_brm2 <- brm(bf(W_G~A*FL_CM^B, A + B ~ 1, nl = TRUE),data = ysdf_2506_FL_W,
prior = ys2506_FLW_prior2,
warmup = 1000,
iter = 5000,
chains = 4,
cores = 4,
control = list(adapt_delta = 0.8))
print(summary(ys2506_FLW_brm2), digits = 6)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: W_G ~ A * FL_CM^B
## A ~ 1
## B ~ 1
## Data: ysdf_2506_FL_W (Number of observations: 84)
## Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 1;
## total post-warmup draws = 16000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## A_Intercept 0.013747 0.002079 0.010006 0.018250 1.000664 3399 2883
## B_Intercept 3.052279 0.041702 2.970542 3.136934 1.000821 3373 2853
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 41.197746 3.279276 35.279404 48.071526 1.000810 6030 5469
##
## 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).Plot a histogram of the residuals to check normality of residuals
hist(resid(ys2506_FLW_brm2))
Calculate R-squared
rss_fitted <- sum(residuals(ys2506_FLW_brm2)^2)
rss_constant <- sum((ysdf_2506_FL_W$W_G - mean(ysdf_2506_FL_W$W_G))^2)
approx_r_squared <- 1 - rss_fitted / rss_constant
approx_r_squared ## [1] 0.829233Calculate the predicted W_G values based on this NLS model
ys2506_FLW_brm2_x <- c(17:46) #rounded actual min and max TL values
ys2506_FLW_brm2_predic <- data.frame(ys2506_FLW_brm2_x)
ys2506_FLW_brm2_f <- function(ys2506_FLW_brm2_x) {
#formula using calculated parameters
return(0.013746*ys2506_FLW_brm2_x^3.052258)
}
ys2506_FLW_brm2_predic$ys2506_FLW_brm2_y <- sapply(ys2506_FLW_brm2_predic$ys2506_FLW_brm2_x, ys2506_FLW_brm2_f)
View(ys2506_FLW_brm2_predic)Plot
ys2506_FLW_brm2_plottemp <- ys2506_FLW_brm1_plot +
geom_smooth(aes(x = ys2506_FLW_brm2_x, y = ys2506_FLW_brm2_y),
data = ys2506_FLW_brm2_predic,
method = "loess",
se = F,
color = safe_colorblind_palette[5],
linewidth = 0.75)
# add label with equation
ys2506_FLW_brm2_plot <- ys2506_FLW_brm2_plottemp +
geom_text(x = 18, y = 1200,
label = expression(W~'(g)' == 9.855%*%10^{-03} * FL~'(cm)'^{3.074}~"," ~ R^{2} == 0.829 ~ '(BRMS)'),
colour = safe_colorblind_palette[5],
hjust = 0)
ys2506_FLW_brm2_plot
##Final decision: NLS
Because the r-squared is larger for the NLS model and the assumptions are met, I will use this model. It is a good sign, however, that the NLS and BRMS parameter estimations are very similar!
Plot the model
ys2506_FLW_nls_plot_title <- expression(italic("L. vivanus")~"FL-W Relationship (n = 84)")
ys2506_FLW_nls_plottemp <- ggplot(ysdf_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 = ys2506_FLW_nls_plot_title) +
theme_classic() +
geom_smooth(aes(x = ys2506_FLW_nls_x, y = ys2506_FLW_nls_y),
data = ys2506_FLW_nls_predic,
method = "loess",
se = F,
color = turbo_colors_onboard[6],
linewidth = 0.75,
alpha = 0.75) +
ylim(0, 1600) +
xlim(12, 48)
# add label with equation
ys2506_FLW_nls_plot <- ys2506_FLW_nls_plottemp +
geom_label(x = 12, y = 1500,
label = expression(W~'(g)' == 1.332%*%10^{-02} * FL~'(cm)'^{3.058}~"," ~ R^{2} == 0.986 ~ '(NLS)'),
fill = "#edde94",
hjust = 0)
ys2506_FLW_nls_plot
Save plot
ggsave("ys2506_FLW_nls_plot.png", plot = ys2506_FLW_nls_plot, width = 7.2, height = 4, dpi = 400)## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_smooth()`).## 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
ysdf_2506_TL_W <- read.csv("ys2506df.csv")
ysdf_2506_TL_W$catch_month <- factor(ysdf_2506_TL_W$catch_month)
ysdf_2506_TL_W$year <- factor(ysdf_2506_TL_W$year)
ysdf_2506_TL_W$species <- factor(ysdf_2506_TL_W$species)
ysdf_2506_TL_W$common_name <- factor(ysdf_2506_TL_W$common_name)
ysdf_2506_TL_W$frozen <- factor(ysdf_2506_TL_W$frozen)
ysdf_2506_TL_W$catch_method <- factor(ysdf_2506_TL_W$catch_method)
ysdf_2506_TL_W$bait <- factor(ysdf_2506_TL_W$bait)
ysdf_2506_TL_W$area <- factor(ysdf_2506_TL_W$area)
ysdf_2506_TL_W$gutted <- factor(ysdf_2506_TL_W$gutted)
ysdf_2506_TL_W$sex <- factor(ysdf_2506_TL_W$sex)
ysdf_2506_TL_W$maturity <- factor(ysdf_2506_TL_W$maturity)
# Remove gutted samples
ysdf_2506_TL_W <- subset(filter(ysdf_2506_TL_W, gutted == "N"))
# remove NA rows
ysdf_2506_TL_W <- ysdf_2506_TL_W%>%
drop_na(TL_CM)
View(ysdf_2506_TL_W)
summary(ysdf_2506_TL_W)## X catch_date catch_month landing_time
## Min. : 1.00 Length:86 Sep :32 Length:86
## 1st Qu.:22.25 Class :character Jun :19 Class :character
## Median :43.50 Mode :character May :14 Mode :character
## Mean :43.50 Jul : 9
## 3rd Qu.:64.75 Mar : 5
## Max. :86.00 : 3
## (Other): 4
## dissection_date dissection_time year month
## Length:86 Length:86 2023:57 Length:86
## Class :character Class :character 2024:29 Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
## observer species common_name
## Length:86 Lutjanus vivanus:86 Yelloweye Snapper:86
## Class :character
## Mode :character
##
##
##
##
## fish_code frozen catch_method bait
## Length:86 N:26 :48 :54
## Class :character Y:60 Dropline : 8 Cowskin : 1
## Mode :character Lobster trap: 1 Japanese bait :23
## Redfish trap:29 Japanese bait; Squid: 8
##
##
##
## area SL_CM FL_CM TL_CM HH_CM
## :25 Min. :16.00 Min. :17.20 Min. :18.30 Min. :3.600
## A3/A4: 5 1st Qu.:24.68 1st Qu.:21.73 1st Qu.:23.20 1st Qu.:5.250
## B4 :13 Median :29.20 Median :26.15 Median :28.60 Median :6.600
## C2 :24 Mean :28.70 Mean :27.57 Mean :29.77 Mean :6.317
## D3 : 6 3rd Qu.:32.10 3rd Qu.:33.08 3rd Qu.:36.08 3rd Qu.:7.450
## D5 :13 Max. :40.40 Max. :45.90 Max. :49.20 Max. :9.400
## NA's :48 NA's :63
## ED_CM W_G finclip otoliths gutted
## Min. :1.20 Min. : 83.9 Length:86 Length:86 N:86
## 1st Qu.:1.65 1st Qu.: 162.7 Class :character Class :character
## Median :1.90 Median : 281.9 Mode :character Mode :character
## Mean :1.87 Mean : 429.0
## 3rd Qu.:2.10 3rd Qu.: 583.3
## Max. :3.00 Max. :1714.7
## NA's :63
## gonads_present sex maturity G1L_CM G2L_CM
## Length:86 :43 :57 Min. : 1.000 Min. : 0.000
## Class :character F:22 A:13 1st Qu.: 2.250 1st Qu.: 2.200
## Mode :character M:21 B:12 Median : 3.800 Median : 3.600
## C: 2 Mean : 4.441 Mean : 4.261
## D: 2 3rd Qu.: 5.950 3rd Qu.: 5.500
## Max. :11.500 Max. :11.600
## NA's :27 NA's :29
## GW_G stomach._sample
## Min. :0.0000 Length:86
## 1st Qu.:0.0000 Class :character
## Median :0.4000 Mode :character
## Mean :0.9695
## 3rd Qu.:1.2500
## Max. :5.6000
## NA's :27Check for outliers in data: calculate z-scores
yszscore_TL_CM <- (ysdf_2506_TL_W$TL_CM - mean(ysdf_2506_TL_W$TL_CM)) / sd(ysdf_2506_TL_W$TL_CM)
ysdf_2506_TL_W$zscore_TL_CM <- yszscore_TL_CM
yszscore_Wgrams <- (ysdf_2506_TL_W$W_G - mean(ysdf_2506_TL_W$W_G)) / sd(ysdf_2506_TL_W$W_G)
ysdf_2506_TL_W$zscore_Wgrams <- yszscore_Wgrams
View(ysdf_2506_TL_W)z-scores above 3 are generally considered to be outliers -> remove from the dataset
Filter the data
ysdf_2506_TL_W <- ysdf_2506_TL_W[ysdf_2506_TL_W$zscore_Wgrams < 3, ]
View(ysdf_2506_TL_W)
summary(ysdf_2506_TL_W)## X catch_date catch_month landing_time
## Min. : 1.00 Length:84 Sep :32 Length:84
## 1st Qu.:23.75 Class :character Jun :18 Class :character
## Median :44.50 Mode :character May :13 Mode :character
## Mean :44.13 Jul : 9
## 3rd Qu.:65.25 Mar : 5
## Max. :86.00 : 3
## (Other): 4
## dissection_date dissection_time year month
## Length:84 Length:84 2023:57 Length:84
## Class :character Class :character 2024:27 Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
## observer species common_name
## Length:84 Lutjanus vivanus:84 Yelloweye Snapper:84
## Class :character
## Mode :character
##
##
##
##
## fish_code frozen catch_method bait
## Length:84 N:24 :48 :54
## Class :character Y:60 Dropline : 8 Cowskin : 1
## Mode :character Lobster trap: 1 Japanese bait :21
## Redfish trap:27 Japanese bait; Squid: 8
##
##
##
## area SL_CM FL_CM TL_CM HH_CM
## :25 Min. :16.00 Min. :17.20 Min. :18.30 Min. :3.600
## A3/A4: 5 1st Qu.:24.32 1st Qu.:21.65 1st Qu.:23.18 1st Qu.:4.900
## B4 :13 Median :28.70 Median :25.55 Median :27.95 Median :6.500
## C2 :24 Mean :28.24 Mean :27.23 Mean :29.51 Mean :6.124
## D3 : 6 3rd Qu.:31.65 3rd Qu.:32.62 3rd Qu.:35.55 3rd Qu.:6.700
## D5 :11 Max. :40.40 Max. :45.90 Max. :49.20 Max. :9.400
## NA's :48 NA's :63
## ED_CM W_G finclip otoliths gutted
## Min. :1.20 Min. : 83.9 Length:84 Length:84 N:84
## 1st Qu.:1.60 1st Qu.: 161.3 Class :character Class :character
## Median :1.90 Median : 272.4 Mode :character Mode :character
## Mean :1.81 Mean : 399.6
## 3rd Qu.:2.10 3rd Qu.: 562.1
## Max. :2.50 Max. :1553.3
## NA's :63
## gonads_present sex maturity G1L_CM G2L_CM
## Length:84 :43 :57 Min. : 1.000 Min. : 0.000
## Class :character F:21 A:12 1st Qu.: 2.200 1st Qu.: 2.050
## Mode :character M:20 B:11 Median : 3.700 Median : 3.200
## C: 2 Mean : 4.195 Mean : 4.002
## D: 2 3rd Qu.: 5.700 3rd Qu.: 5.300
## Max. :10.900 Max. :11.100
## NA's :27 NA's :29
## GW_G stomach._sample zscore_TL_CM zscore_Wgrams
## Min. :0.0000 Length:84 Min. :-1.44840 Min. :-0.88288
## 1st Qu.:0.0000 Class :character 1st Qu.:-0.83254 1st Qu.:-0.68482
## Median :0.3000 Mode :character Median :-0.22931 Median :-0.40054
## Mean :0.8807 Mean :-0.03214 Mean :-0.07534
## 3rd Qu.:1.1000 3rd Qu.: 0.73081 3rd Qu.: 0.34036
## Max. :5.6000 Max. : 2.45523 Max. : 2.87605
## NA's :27First look at the data with a scatterplot
fill = sex
ggplot(ysdf_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 vivanus, n = 84") +
scale_fill_viridis(option = "D", discrete = TRUE)
Linear model
ys2506_TLW_lm1 <- lm(TL_CM ~ W_G, data = ysdf_2506_TL_W)
summary(ys2506_TLW_lm1)##
## Call:
## lm(formula = TL_CM ~ W_G, data = ysdf_2506_TL_W)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.952 -1.646 0.225 1.656 4.149
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.084e+01 3.897e-01 53.47 <2e-16 ***
## W_G 2.170e-02 7.405e-04 29.31 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.324 on 82 degrees of freedom
## Multiple R-squared: 0.9129, Adjusted R-squared: 0.9118
## F-statistic: 859.1 on 1 and 82 DF, p-value: < 2.2e-16Scatterplot the residuals against the fitted values to check homogeneity of variance
plot(resid(ys2506_TLW_lm1) ~ fitted(ys2506_TLW_lm1))
looks nonlinear
Plot a histogram to check normality of residuals
qqnorm(resid(ys2506_TLW_lm1))
qqline(resid(ys2506_TLW_lm1))
Plot the model
ys2506_TLW_lm1_plottemp <- ggplot(ysdf_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 = "Lutjanis vivanus, n = 84") +
theme_classic() +
geom_smooth(method = lm,
aes(TL_CM, W_G),
linewidth = 0.5,
se = FALSE,
color = safe_colorblind_palette[1])
# add label with equation
ys2506_TLW_lm1_plot <- ys2506_TLW_lm1_plottemp +
geom_text(x = 18, y = 1500,
label = expression(W~'(g)' == 0.022 * TL~'(cm)' + 20.84~"," ~ r^{2} == 0.9118 ~ "(Linear)"),
colour = safe_colorblind_palette[1],
hjust = 0)
ys2506_TLW_lm1_plot
NLS
ys2506_TLW_nls <- nls_multstart(W_G~A*TL_CM^B,
data = ysdf_2506_TL_W,
iter = 500,
start_lower = c(A=0, B=2),
start_upper = c(A=0.5, B=3.5))
summary(ys2506_TLW_nls)##
## Formula: W_G ~ A * TL_CM^B
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## A 0.009313 0.001458 6.386 9.67e-09 ***
## B 3.086674 0.042263 73.035 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 38 on 82 degrees of freedom
##
## Number of iterations to convergence: 10
## Achieved convergence tolerance: 1.49e-08Plot a histogram of the residuals to check normality of residuals
qqnorm(resid(ys2506_TLW_nls))
qqline(resid(ys2506_TLW_nls))
Calculate R-squared
rss_fitted <- sum(residuals(ys2506_TLW_nls)^2)
rss_constant <- sum((ysdf_2506_TL_W$W_G - mean(ysdf_2506_TL_W$W_G))^2)
approx_r_squared <- 1 - rss_fitted / rss_constant
approx_r_squared## [1] 0.9879802Calculate AIC
AIC(ys2506_TLW_nls)## [1] 853.4719Calculate the predicted W_kg values based on this NLS model
ys2506_TLW_nls_x <- c(18:50) #rounded actual min and max TL_MM values
ys2506_TLW_nls_predic <- data.frame(ys2506_TLW_nls_x)
ys2506_TLW_nls_f <- function(ys2506_TLW_nls_x) {
#formula using calculated parameters
return(0.009313*ys2506_TLW_nls_x^3.086674)
}
ys2506_TLW_nls_predic$ys2506_TLW_nls_y <- sapply(ys2506_TLW_nls_predic$ys2506_TLW_nls_x, ys2506_TLW_nls_f)
View(ys2506_TLW_nls_predic)Plot the model on the same graph as the other model
ys2506_TLW_nls_plottemp <- ys2506_TLW_lm1_plot +
geom_smooth(aes(x = ys2506_TLW_nls_x, y = ys2506_TLW_nls_y),
data = ys2506_TLW_nls_predic,
method = "loess",
se = F,
color = safe_colorblind_palette[2],
linewidth = 0.75)
# add label with equation
ys2506_TLW_nls_plot <- ys2506_TLW_nls_plottemp +
geom_text(x = 18, y = 1400,
label = expression(W~'(g)' == 9.313%*%10^{-03} * TL~'(cm)'^{3.087}~"," ~ R^{2} == 0.985 ~ '(NLS)'),
colour = safe_colorblind_palette[2],
hjust = 0)
ys2506_TLW_nls_plot
Brms: with nls parameters to start
set priors and run model
#set priors
ys2506_TLW_prior1 <- prior(normal(0.009313, 0.001458) ,nlpar = "A")+prior(normal(3.086674, 0.042263),nlpar = "B")
ys2506_TLW_brm1 <- brm(bf(W_G~A*TL_CM^B, A + B ~ 1, nl = TRUE),data = ysdf_2506_TL_W,
prior = ys2506_TLW_prior1,
warmup = 1000,
iter = 5000,
chains = 4,
cores = 4,
control = list(adapt_delta = 0.8))
print(summary(ys2506_TLW_brm1), digits = 6)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: W_G ~ A * TL_CM^B
## A ~ 1
## B ~ 1
## Data: ysdf_2506_TL_W (Number of observations: 84)
## Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 1;
## total post-warmup draws = 16000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## A_Intercept 0.009400 0.000833 0.007845 0.011082 1.000731 5288 5454
## B_Intercept 3.085172 0.023980 3.039996 3.133313 1.000748 5277 5099
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 38.420445 3.037137 32.887421 44.863338 1.000467 5524 5614
##
## 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).Plot a histogram of the residuals to check normality of residuals
hist(resid(ys2506_TLW_brm1))
Calculate R-squared
rss_fitted <- sum(residuals(ys2506_TLW_brm1)^2)
rss_constant <- sum((ysdf_2506_TL_W$W_G - mean(ysdf_2506_TL_W$W_G))^2)
approx_r_squared <- 1 - rss_fitted / rss_constant
approx_r_squared ## [1] 0.852279Calculate the predicted W_G values based on this NLS model
ys2506_TLW_brm1_x <- c(18:50) #rounded actual min and max TL values
ys2506_TLW_brm1_predic <- data.frame(ys2506_TLW_brm1_x)
ys2506_TLW_brm1_f <- function(ys2506_TLW_brm1_x) {
#formula using calculated parameters
return(0.009407*ys2506_TLW_brm1_x^3.085021)
}
ys2506_TLW_brm1_predic$ys2506_TLW_brm1_y <- sapply(ys2506_TLW_brm1_predic$ys2506_TLW_brm1_x, ys2506_TLW_brm1_f)
View(ys2506_TLW_brm1_predic)Plot
ys2506_TLW_brm1_plottemp <- ys2506_TLW_nls_plot +
geom_smooth(aes(x = ys2506_TLW_brm1_x, y = ys2506_TLW_brm1_y),
data = ys2506_TLW_brm1_predic,
method = "loess",
se = F,
color = safe_colorblind_palette[4],
linewidth = 0.75)
# add label with equation
ys2506_TLW_brm1_plot <- ys2506_TLW_brm1_plottemp +
geom_text(x = 18, y = 1300,
label = expression(W~'(g)' == 9.407%*%10^{-03} * TL~'(cm)'^{3.085}~"," ~ R^{2} == 0.852 ~ '(BRMS)'),
colour = safe_colorblind_palette[4],
hjust = 0)
ys2506_TLW_brm1_plot
The NLS model and the BRMS models are very similar which is a good sign
Use the NLS model since the R-squared is better?
Can try getting a brms model to fit even better
Brms2: with previous brms parameters as priors
also tightening the normal distribution around the starting priors using the difference between the upper and lower CI from the previous models
#A UCI-LCI
0.011192-0.007815 ## [1] 0.003377#B UCI-LCI
3.133829-3.036511 ## [1] 0.097318set priors and run model
#set priors
ys2506_TLW_prior2 <- prior(normal(0.009407, 0.003377) ,nlpar = "A")+prior(normal(3.085021, 0.097318),nlpar = "B")
ys2506_TLW_brm2 <- brm(bf(W_G~A*TL_CM^B, A + B ~ 1, nl = TRUE), data = ysdf_2506_TL_W,
prior = ys2506_TLW_prior2,
warmup = 1000,
iter = 5000,
chains = 4,
cores = 4,
control = list(adapt_delta = 0.8))
print(summary(ys2506_TLW_brm2), digits = 6)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: W_G ~ A * TL_CM^B
## A ~ 1
## B ~ 1
## Data: ysdf_2506_TL_W (Number of observations: 84)
## Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 1;
## total post-warmup draws = 16000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## A_Intercept 0.009556 0.001268 0.007236 0.012250 1.000849 5180 5135
## B_Intercept 3.082080 0.035901 3.012800 3.154857 1.000850 5165 5002
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 38.508926 3.055624 33.051328 44.995125 1.000865 5259 5322
##
## 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).Plot a histogram of the residuals to check normality of residuals
hist(resid(ys2506_TLW_brm2))
Calculate R-squared
rss_fitted <- sum(residuals(ys2506_TLW_brm2)^2)
rss_constant <- sum((ysdf_2506_TL_W$W_G - mean(ysdf_2506_TL_W$W_G))^2)
approx_r_squared <- 1 - rss_fitted / rss_constant
approx_r_squared ## [1] 0.8505885Calculate the predicted W_G values based on this NLS model
ys2506_TLW_brm2_x <- c(18:50) #rounded actual min and max TL values
ys2506_TLW_brm2_predic <- data.frame(ys2506_TLW_brm2_x)
ys2506_TLW_brm2_f <- function(ys2506_TLW_brm2_x) {
#formula using calculated parameters
return(0.009543*ys2506_TLW_brm2_x^3.082496)
}
ys2506_TLW_brm2_predic$ys2506_TLW_brm2_y <- sapply(ys2506_TLW_brm2_predic$ys2506_TLW_brm2_x, ys2506_TLW_brm2_f)
View(ys2506_TLW_brm2_predic)Plot
ys2506_TLW_brm2_plottemp <- ys2506_TLW_brm1_plot +
geom_smooth(aes(x = ys2506_TLW_brm2_x, y = ys2506_TLW_brm2_y),
data = ys2506_TLW_brm2_predic,
method = "loess",
se = F,
color = safe_colorblind_palette[5],
linewidth = 0.75,
linetype = 2)
# add label with equation
ys2506_TLW_brm2_plot <- ys2506_TLW_brm2_plottemp +
geom_text(x = 18, y = 1200,
label = expression(W~'(g)' == 9.543%*%10^{-03} * TL~'(cm)'^{3.082}~"," ~ R^{2} == "0.850" ~ '(BRMS)'),
colour = safe_colorblind_palette[5],
hjust = 0)
ys2506_TLW_brm2_plot
##Final decision: NLS
the second brms has a better R-squared than the first but the NLS is still better. It is a good sign that the parameters are quite similar in all three nonlinear regression growth functions.
Load data file
ysdf_2506_SL_W <- read.csv("ys2506df.csv")
ysdf_2506_SL_W$catch_month <- factor(ysdf_2506_SL_W$catch_month)
ysdf_2506_SL_W$year <- factor(ysdf_2506_SL_W$year)
ysdf_2506_SL_W$species <- factor(ysdf_2506_SL_W$species)
ysdf_2506_SL_W$common_name <- factor(ysdf_2506_SL_W$common_name)
ysdf_2506_SL_W$frozen <- factor(ysdf_2506_SL_W$frozen)
ysdf_2506_SL_W$catch_method <- factor(ysdf_2506_SL_W$catch_method)
ysdf_2506_SL_W$bait <- factor(ysdf_2506_SL_W$bait)
ysdf_2506_SL_W$area <- factor(ysdf_2506_SL_W$area)
ysdf_2506_SL_W$gutted <- factor(ysdf_2506_SL_W$gutted)
ysdf_2506_SL_W$sex <- factor(ysdf_2506_SL_W$sex)
ysdf_2506_SL_W$maturity <- factor(ysdf_2506_SL_W$maturity)
# Remove gutted samples
ysdf_2506_SL_W <- subset(filter(ysdf_2506_SL_W, gutted == "N"))
# remove NA rows
ysdf_2506_SL_W <- ysdf_2506_SL_W%>%
drop_na(SL_CM)
View(ysdf_2506_SL_W)
summary(ysdf_2506_SL_W)## X catch_date catch_month landing_time
## Min. : 1.00 Length:38 May :14 Length:38
## 1st Qu.:10.25 Class :character Jun : 9 Class :character
## Median :19.50 Mode :character Sep : 8 Mode :character
## Mean :19.50 Mar : 5
## 3rd Qu.:28.75 Apr : 1
## Max. :38.00 Dec : 1
## (Other): 0
## dissection_date dissection_time year month
## Length:38 Length:38 2023: 9 Length:38
## Class :character Class :character 2024:29 Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
## observer species common_name
## Length:38 Lutjanus vivanus:38 Yelloweye Snapper:38
## Class :character
## Mode :character
##
##
##
##
## fish_code frozen catch_method bait
## Length:38 N:26 : 0 : 6
## Class :character Y:12 Dropline : 8 Cowskin : 1
## Mode :character Lobster trap: 1 Japanese bait :23
## Redfish trap:29 Japanese bait; Squid: 8
##
##
##
## area SL_CM FL_CM TL_CM HH_CM
## : 1 Min. :16.00 Min. :18.90 Min. :20.20 Min. :3.600
## A3/A4: 5 1st Qu.:24.68 1st Qu.:28.25 1st Qu.:30.07 1st Qu.:5.250
## B4 :13 Median :29.20 Median :33.45 Median :36.40 Median :6.600
## C2 : 0 Mean :28.70 Mean :32.89 Mean :35.38 Mean :6.317
## D3 : 6 3rd Qu.:32.10 3rd Qu.:36.73 3rd Qu.:39.95 3rd Qu.:7.450
## D5 :13 Max. :40.40 Max. :45.90 Max. :49.20 Max. :9.400
## NA's :15
## ED_CM W_G finclip otoliths gutted
## Min. :1.20 Min. : 106.8 Length:38 Length:38 N:38
## 1st Qu.:1.65 1st Qu.: 371.1 Class :character Class :character
## Median :1.90 Median : 621.8 Mode :character Mode :character
## Mean :1.87 Mean : 693.8
## 3rd Qu.:2.10 3rd Qu.: 855.2
## Max. :3.00 Max. :1714.7
## NA's :15
## gonads_present sex maturity G1L_CM G2L_CM
## Length:38 : 2 : 9 Min. : 1.000 Min. : 1.000
## Class :character F:21 A:13 1st Qu.: 2.925 1st Qu.: 2.700
## Mode :character M:15 B:12 Median : 4.200 Median : 4.000
## C: 2 Mean : 5.129 Mean : 5.186
## D: 2 3rd Qu.: 6.875 3rd Qu.: 7.000
## Max. :11.500 Max. :11.600
## NA's :1
## GW_G stomach._sample
## Min. :0.000 Length:38
## 1st Qu.:0.200 Class :character
## Median :0.950 Mode :character
## Mean :1.397
## 3rd Qu.:1.700
## Max. :5.600
## Check for outliers in data: calculate z-scores
yszscore_SL_CM <- (ysdf_2506_SL_W$SL_CM - mean(ysdf_2506_SL_W$SL_CM)) / sd(ysdf_2506_SL_W$SL_CM)
ysdf_2506_SL_W$zscore_SL_CM <- yszscore_SL_CM
yszscore_Wgrams <- (ysdf_2506_SL_W$W_G - mean(ysdf_2506_SL_W$W_G)) / sd(ysdf_2506_SL_W$W_G)
ysdf_2506_SL_W$zscore_Wgrams <- yszscore_Wgrams
View(ysdf_2506_SL_W)No z-scores above 3
First look at the data with a scatterplot
fill = sex
ggplot(ysdf_2506_SL_W, aes(x = SL_CM, y = W_G)) +
geom_point(aes(fill = sex), color = "#1B1B1B", pch = 21, size = 2) +
theme_classic() +
labs(x = "Standard Length in cm", y = "Weight in grams", title = "Lutjanus vivanus, n = 38") +
scale_fill_viridis(option = "D", discrete = TRUE)
LM Linear model FL-W
ys2506_FLW_lm <- lm(SL_CM ~ W_G, data = ysdf_2506_SL_W)
summary(ys2506_FLW_lm)##
## Call:
## lm(formula = SL_CM ~ W_G, data = ysdf_2506_SL_W)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.1986 -1.1175 0.9543 1.5461 2.3672
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.946e+01 6.759e-01 28.78 <2e-16 ***
## W_G 1.332e-02 8.188e-04 16.27 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.258 on 36 degrees of freedom
## Multiple R-squared: 0.8803, Adjusted R-squared: 0.877
## F-statistic: 264.7 on 1 and 36 DF, p-value: < 2.2e-16Scatterplot the residuals against the fitted values to check homogeneity of variance
plot(resid(ys2506_FLW_lm) ~ fitted(ys2506_FLW_lm))
looks like it indicates a nonlinear relationship
Plot a histogram to check normality of residuals
qqnorm(resid(ys2506_FLW_lm))
qqline(resid(ys2506_FLW_lm))
not normal
Plot the model
ys2506_SLW_lm_plottemp <- ggplot(ysdf_2506_SL_W, aes(x = SL_CM, y = W_G)) +
geom_point(aes(SL_CM, W_G),
color = "black",
alpha = 0.6) +
labs(x = "Standard Length in cm",
y = "Weight in grams",
title = "Lutjanus vivanus, n = 38") +
theme_classic() +
geom_smooth(method = lm,
aes(SL_CM, W_G),
linewidth = 0.5,
se = FALSE,
color = safe_colorblind_palette[1])
# add label with equation
ys2506_SLW_lm_plot <- ys2506_SLW_lm_plottemp +
geom_text(x = 16, y = 1500,
label = expression(W~'(g)' == 0.029 * SL~'(cm)' + 13.611~"," ~ r^{2} == 0.9529 ~ "(Linear)"),
colour = safe_colorblind_palette[1],
hjust = 0)
ys2506_SLW_lm_plot
NLS Models
NLS 1
ys2506_SLW_nls <- nls_multstart(W_G~A*SL_CM^B,
data = ysdf_2506_SL_W,
iter = 500,
start_lower = c(A=0, B=2),
start_upper = c(A=0.5, B=3.5))
print(summary(ys2506_SLW_nls), digits = 4)##
## Formula: W_G ~ A * SL_CM^B
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## A 0.01869 0.01093 1.71 0.0959 .
## B 3.09230 0.16492 18.75 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 116.6 on 36 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
#check residuals
qqnorm(resid(ys2506_SLW_nls))
qqline(resid(ys2506_SLW_nls))
okay
Calculate R-squared
rss_fitted <- sum(residuals(ys2506_SLW_nls)^2)
rss_constant <- sum((ysdf_2506_SL_W$W_G - mean(ysdf_2506_SL_W$W_G))^2)
approx_r_squared <- 1 - rss_fitted / rss_constant
approx_r_squared## [1] 0.9357097Calculate AIC
AIC(ys2506_SLW_nls)## [1] 473.4177Calculate the predicted W_g values based on this NLS model
ys2506_SLW_nls_x <- c(16:41) #rounded actual min and max FL_MM values
ys2506_SLW_nls_predic <- data.frame(ys2506_SLW_nls_x)
ys2506_SLW_nls_f <- function(ys2506_SLW_nls_x) {
#formula using calculated parameters
return(0.01869*ys2506_SLW_nls_x^3.09230)
}
ys2506_SLW_nls_predic$ys2506_SLW_nls_y <- sapply(ys2506_SLW_nls_predic$ys2506_SLW_nls_x, ys2506_SLW_nls_f)
View(ys2506_SLW_nls_predic)Plot the model on the same graph as the other model
ys2506_SLW_nls_plottemp <- ys2506_SLW_lm_plot +
geom_smooth(aes(x = ys2506_SLW_nls_x, y = ys2506_SLW_nls_y),
data = ys2506_SLW_nls_predic,
method = "loess",
se = F,
color = safe_colorblind_palette[2],
linewidth = 0.75,
alpha = 0.75)
# add label with equation
ys2506_SLW_nls_plot <- ys2506_SLW_nls_plottemp +
geom_text(x = 16, y = 1350,
label = expression(W~'(g)' == 1.869%*%10^{-02} * SL~'(cm)'^{3.092}~"," ~ R^{2} == 0.979 ~ '(NLS)'),
colour = safe_colorblind_palette[2],
hjust = 0)
ys2506_SLW_nls_plot
NLS 2
#A upper = estimate + SE
0.01869+0.01093## [1] 0.02962#A lower = estimate - SE
0.01869-0.01093## [1] 0.00776#B upper = estimate + SE
3.09230+0.16492## [1] 3.25722#B lower = estimate - SE
3.09230-0.16492## [1] 2.92738ys2506_SLW_nls2 <- nls_multstart(W_G~A*SL_CM^B,
data = ysdf_2506_SL_W,
iter = 500,
start_lower = c(A=0.00776, B=2.92738),
start_upper = c(A=0.02962, B=3.25722))
print(summary(ys2506_SLW_nls2), digits = 4)##
## Formula: W_G ~ A * SL_CM^B
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## A 0.01869 0.01093 1.71 0.0959 .
## B 3.09230 0.16492 18.75 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 116.6 on 36 degrees of freedom
##
## Number of iterations to convergence: 9
## Achieved convergence tolerance: 1.49e-08parameter estimates are the same as before with same SE so that is a good sign
Brms: with nls parameters to start
set priors and run model
#set priors
ys2506_SLW_prior1 <- prior(normal(0.01869, 0.01093) ,nlpar = "A")+prior(normal(3.09230, 0.16492),nlpar = "B")
ys2506_SLW_brm1 <- brm(bf(W_G~A*SL_CM^B, A + B ~ 1, nl = TRUE),data = ysdf_2506_SL_W,
prior = ys2506_SLW_prior1,
warmup = 1000,
iter = 5000,
chains = 4,
cores = 4,
control = list(adapt_delta = 0.8))
print(summary(ys2506_SLW_brm1), digits = 6)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: W_G ~ A * SL_CM^B
## A ~ 1
## B ~ 1
## Data: ysdf_2506_SL_W (Number of observations: 38)
## Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 1;
## total post-warmup draws = 16000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## A_Intercept 0.020609 0.006298 0.010087 0.034505 1.000336 3302 1919
## B_Intercept 3.077746 0.088250 2.919573 3.266601 1.000352 3319 1820
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 119.238312 14.393477 94.966608 151.720823 1.000880 4156 4885
##
## 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).Plot a histogram of the residuals to check normality of residuals
hist(resid(ys2506_SLW_brm1))
Calculate R-squared
rss_fitted <- sum(residuals(ys2506_SLW_brm1)^2)
rss_constant <- sum((ysdf_2506_SL_W$W_G - mean(ysdf_2506_SL_W$W_G))^2)
approx_r_squared <- 1 - rss_fitted / rss_constant
approx_r_squared ## [1] 0.1559244Calculate the predicted W_G values based on this NLS model
ys2506_SLW_brm1_x <- c(16:41) #rounded actual min and max TL values
ys2506_SLW_brm1_predic <- data.frame(ys2506_SLW_brm1_x)
ys2506_SLW_brm1_f <- function(ys2506_SLW_brm1_x) {
#formula using calculated parameters
return(0.020679*ys2506_SLW_brm1_x^3.076635)
}
ys2506_SLW_brm1_predic$ys2506_SLW_brm1_y <- sapply(ys2506_SLW_brm1_predic$ys2506_SLW_brm1_x, ys2506_SLW_brm1_f)
View(ys2506_SLW_brm1_predic)Plot
ys2506_SLW_brm1_plottemp <- ys2506_SLW_nls_plot +
geom_smooth(aes(x = ys2506_SLW_brm1_x, y = ys2506_SLW_brm1_y),
data = ys2506_SLW_brm1_predic,
method = "loess",
se = F,
color = safe_colorblind_palette[3],
linewidth = 0.75)
# add label with equation
ys2506_SLW_brm1_plot <- ys2506_SLW_brm1_plottemp +
geom_text(x = 16, y = 1200,
label = expression(W~'(g)' == 2.068%*%10^{-02} * SL~'(cm)'^{3.077}~"," ~ R^{2} == 0.154 ~ '(BRMS)'),
colour = safe_colorblind_palette[3],
hjust = 0)
ys2506_SLW_brm1_plot
##Final decision: NLS
Plot the model
ys2506_SLW_nls_plot_title <- expression(italic("L. vivanus") ~ "SL-W Relationship (n = 38)")
ys2506_sqSLW_nls_plottemp <- ggplot(ysdf_2506_SL_W, aes(x = SL_CM, y = W_G)) +
geom_point(aes(SL_CM, W_G),
color = "#242124",
alpha = 0.8) +
labs(x = "Standard Length in cm",
y = "Weight in grams",
title = ys2506_SLW_nls_plot_title) +
theme_classic() +
geom_smooth(aes(x = ys2506_SLW_nls_x, y = ys2506_SLW_nls_y),
data = ys2506_SLW_nls_predic,
method = "loess",
se = F,
color = turbo_colors_onboard[6],
linewidth = 0.75,
alpha = 0.75)
# add label with equation
ys2506_sqSLW_nls_plot <- ys2506_sqSLW_nls_plottemp +
geom_text(x = 16, y = 1800,
label = expression(W~'(g)' == 3.347%*%10^{-02} * SL~'(cm)'^{2.895}~"," ~ R^{2} == 0.979 ~ '(NLS)'),
colour = turbo_colors_onboard[6],
hjust = 0)
ys2506_sqSLW_nls_plot
Maturity
Load data file
ysdf_2506_maturity <- read.csv("ys2506df.csv")
ysdf_2506_maturity$catch_month <- factor(ysdf_2506_maturity$catch_month)
ysdf_2506_maturity$year <- factor(ysdf_2506_maturity$year)
ysdf_2506_maturity$species <- factor(ysdf_2506_maturity$species)
ysdf_2506_maturity$common_name <- factor(ysdf_2506_maturity$common_name)
ysdf_2506_maturity$frozen <- factor(ysdf_2506_maturity$frozen)
ysdf_2506_maturity$catch_method <- factor(ysdf_2506_maturity$catch_method)
ysdf_2506_maturity$bait <- factor(ysdf_2506_maturity$bait)
ysdf_2506_maturity$area <- factor(ysdf_2506_maturity$area)
ysdf_2506_maturity$gutted <- factor(ysdf_2506_maturity$gutted)
ysdf_2506_maturity$sex <- factor(ysdf_2506_maturity$sex)
ysdf_2506_maturity$maturity <- factor(ysdf_2506_maturity$maturity)
#drop the columns with neither sex nor maturity measurements
ysdf_2506_maturity <- ysdf_2506_maturity %>%
filter(maturity != "" | sex != "")calculate sex ratio
ys_totalfemale <- ysdf_2506_maturity %>%
filter(sex == "F") %>%
count()
ys_totalmale <- ysdf_2506_maturity %>%
filter(sex == "M") %>%
count()
ys_sexratio <- ys_totalfemale/ys_totalmale
ys_sexratio #female number## n
## 1 1.047619Descriptive statistics
ys_sex_summary <- ysdf_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)
)
ys_sex_summary## # A tibble: 3 × 8
## sex count mean sd median min max IQR
## <fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 "" 1 18.9 NA 18.9 18.9 18.9 0
## 2 "F" 22 32.5 7.23 32.5 19.2 43.9 9.73
## 3 "M" 21 30.2 8.50 32.7 18.8 45.9 13.6ysdf_2506_maturity <- ysdf_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))
ys_maturity_summary <- ysdf_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)
)
ys_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 13 26.9 7.71 24.6 18.9 43.9 8.7
## 2 mature 16 35.4 5.08 34.2 27.3 44.1 4.65
## 3 <NA> 15 30.2 8.96 31.7 18.8 45.9 14.8