P. aquilonaris body size models
Naomi Bakken
2024-08-15
This R markdown file includes all the R code used to analyze Wenchman Snapper (Pristipomoides aquilonaris) 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
wndf_2506_LLrelationships <- read.csv("wrn2506df.csv")
wndf_2506_LLrelationships$catch_month <- factor(wndf_2506_LLrelationships$catch_month)
wndf_2506_LLrelationships$year <- factor(wndf_2506_LLrelationships$year)
wndf_2506_LLrelationships$species <- factor(wndf_2506_LLrelationships$species)
wndf_2506_LLrelationships$common_name <- factor(wndf_2506_LLrelationships$common_name)
wndf_2506_LLrelationships$frozen <- factor(wndf_2506_LLrelationships$frozen)
wndf_2506_LLrelationships$catch_method <- factor(wndf_2506_LLrelationships$catch_method)
wndf_2506_LLrelationships$bait <- factor(wndf_2506_LLrelationships$bait)
wndf_2506_LLrelationships$area <- factor(wndf_2506_LLrelationships$area)
wndf_2506_LLrelationships$gutted <- factor(wndf_2506_LLrelationships$gutted)
wndf_2506_LLrelationships$sex <- factor(wndf_2506_LLrelationships$sex)
wndf_2506_LLrelationships$maturity <- factor(wndf_2506_LLrelationships$maturity)
summary(wndf_2506_LLrelationships)## X catch_date catch_month landing_time
## Min. : 1.00 Length:24 Apr: 1 Length:24
## 1st Qu.: 6.75 Class :character Jun: 6 Class :character
## Median :12.50 Mode :character May: 1 Mode :character
## Mean :12.50 Sep:16
## 3rd Qu.:18.25
## Max. :24.00
##
## dissection_date dissection_time year month
## Length:24 Length:24 2023:16 Length:24
## Class :character Class :character 2024: 8 Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
## observer species common_name
## Length:24 Pristipomoides aquilonaris:24 Wenchman Snapper:24
## Class :character
## Mode :character
##
##
##
##
## fish_code frozen catch_method bait area
## Length:24 N: 8 Dropline :16 Japanese bait : 8 : 1
## Class :character Y:16 Redfish trap: 8 Japanese bait; Squid:16 B4:16
## Mode :character D4: 1
## D5: 6
##
##
##
## SL_CM FL_CM TL_CM HH_CM ED_CM
## Min. :19.50 Min. :21.60 Min. :25.50 Min. :4.30 Min. :1.500
## 1st Qu.:24.75 1st Qu.:27.27 1st Qu.:32.35 1st Qu.:5.05 1st Qu.:1.700
## Median :28.00 Median :31.15 Median :35.20 Median :5.70 Median :2.000
## Mean :28.10 Mean :31.10 Mean :35.53 Mean :5.60 Mean :2.057
## 3rd Qu.:31.43 3rd Qu.:34.60 3rd Qu.:38.70 3rd Qu.:6.05 3rd Qu.:2.450
## Max. :36.40 Max. :40.30 Max. :45.40 Max. :7.00 Max. :2.600
## NA's :1 NA's :17 NA's :17
## W_G finclip otoliths gutted
## Min. : 197.3 Length:24 Length:24 N:24
## 1st Qu.: 389.6 Class :character Class :character
## Median : 591.7 Mode :character Mode :character
## Mean : 626.0
## 3rd Qu.: 798.0
## Max. :1302.6
##
## gonads_present sex maturity G1L_CM G2L_CM
## Length:24 : 1 :16 Min. :1.500 Min. :2.100
## Class :character F:15 A: 3 1st Qu.:3.175 1st Qu.:3.350
## Mode :character M: 8 B: 5 Median :4.100 Median :3.900
## Mean :3.962 Mean :4.239
## 3rd Qu.:4.625 3rd Qu.:5.000
## Max. :7.000 Max. :7.100
## NA's :1
## GW_G stomach._sample
## Min. : 0.000 Length:24
## 1st Qu.: 0.550 Class :character
## Median : 0.850 Mode :character
## Mean : 3.008
## 3rd Qu.: 2.200
## Max. :35.300
## View(wndf_2506_LLrelationships)
sumtable(wndf_2506_LLrelationships, digits = 4)| Variable | N | Mean | Std. Dev. | Min | Pctl. 25 | Pctl. 75 | Max |
|---|---|---|---|---|---|---|---|
| X | 24 | 12.5 | 7.071 | 1 | 6.75 | 18.25 | 24 |
| catch_date | 24 | ||||||
| … 02/04/2024 | 1 | 4.17% | |||||
| … 05/06/2024 | 5 | 20.83% | |||||
| … 19/06/2024 | 1 | 4.17% | |||||
| … 19/09/2023 | 16 | 66.67% | |||||
| … 29/05/2024 | 1 | 4.17% | |||||
| catch_month | 24 | ||||||
| … Apr | 1 | 4.17% | |||||
| … Jun | 6 | 25% | |||||
| … May | 1 | 4.17% | |||||
| … Sep | 16 | 66.67% | |||||
| landing_time | 24 | ||||||
| … | 16 | 66.67% | |||||
| … 12:45 | 1 | 4.17% | |||||
| … 13:30 | 1 | 4.17% | |||||
| … 14:10 | 1 | 4.17% | |||||
| … 14:45 | 5 | 20.83% | |||||
| year | 24 | ||||||
| … 2023 | 16 | 66.67% | |||||
| … 2024 | 8 | 33.33% | |||||
| month | 24 | ||||||
| … April | 1 | 4.17% | |||||
| … Jun | 6 | 25% | |||||
| … May | 1 | 4.17% | |||||
| … Nov | 16 | 66.67% | |||||
| observer | 24 | ||||||
| … LJ | 16 | 66.67% | |||||
| … NB | 2 | 8.33% | |||||
| … SH | 4 | 16.67% | |||||
| … SM | 2 | 8.33% | |||||
| species | 24 | ||||||
| … Pristipomoides aquilonaris | 24 | 100% | |||||
| common_name | 24 | ||||||
| … Wenchman Snapper | 24 | 100% | |||||
| frozen | 24 | ||||||
| … N | 8 | 33.33% | |||||
| … Y | 16 | 66.67% | |||||
| catch_method | 24 | ||||||
| … Dropline | 16 | 66.67% | |||||
| … Redfish trap | 8 | 33.33% | |||||
| bait | 24 | ||||||
| … Japanese bait | 8 | 33.33% | |||||
| … Japanese bait; Squid | 16 | 66.67% | |||||
| area | 24 | ||||||
| … | 1 | 4.17% | |||||
| … B4 | 16 | 66.67% | |||||
| … D4 | 1 | 4.17% | |||||
| … D5 | 6 | 25% | |||||
| SL_CM | 24 | 28.1 | 4.24 | 19.5 | 24.75 | 31.42 | 36.4 |
| FL_CM | 24 | 31.1 | 4.692 | 21.6 | 27.28 | 34.6 | 40.3 |
| TL_CM | 23 | 35.53 | 4.62 | 25.5 | 32.35 | 38.7 | 45.4 |
| HH_CM | 7 | 5.6 | 0.8981 | 4.3 | 5.05 | 6.05 | 7 |
| ED_CM | 7 | 2.057 | 0.4467 | 1.5 | 1.7 | 2.45 | 2.6 |
| W_G | 24 | 626 | 277.4 | 197.3 | 389.6 | 798 | 1303 |
| finclip | 24 | ||||||
| … Y | 24 | 100% | |||||
| otoliths | 24 | ||||||
| … Y | 24 | 100% | |||||
| gutted | 24 | ||||||
| … N | 24 | 100% | |||||
| gonads_present | 24 | ||||||
| … Y | 24 | 100% | |||||
| sex | 24 | ||||||
| … | 1 | 4.17% | |||||
| … F | 15 | 62.5% | |||||
| … M | 8 | 33.33% | |||||
| maturity | 24 | ||||||
| … | 16 | 66.67% | |||||
| … A | 3 | 12.5% | |||||
| … B | 5 | 20.83% | |||||
| G1L_CM | 24 | 3.962 | 1.336 | 1.5 | 3.175 | 4.625 | 7 |
| G2L_CM | 23 | 4.239 | 1.369 | 2.1 | 3.35 | 5 | 7.1 |
| GW_G | 24 | 3.008 | 7.152 | 0 | 0.55 | 2.2 | 35.3 |
| stomach._sample | 24 | ||||||
| … N | 6 | 25% | |||||
| … Y | 18 | 75% |
Fork length (FL_CM) - Total length (TL_CM) relationship
wn2506_FLTL_plotbase <- ggplot(wndf_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 = "Pristipomoides aquilonaris FL-TL relationship (n = 23)",
fill = "Sex") +
scale_fill_viridis(option = "C",
discrete = TRUE,
begin = 0.2,
end = 0.8,
labels = c("N/A", "Female", "Male"))
wn2506_FLTL_plotbase
Linear model
wn2506_FLTL_lm <- lm(FL_CM ~ TL_CM, data = wndf_2506_LLrelationships)
summary(wn2506_FLTL_lm)##
## Call:
## lm(formula = FL_CM ~ TL_CM, data = wndf_2506_LLrelationships)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.86051 -0.30147 0.00343 0.19328 1.01474
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.62121 0.73694 -2.20 0.0391 *
## TL_CM 0.93232 0.02057 45.32 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4458 on 21 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.9899, Adjusted R-squared: 0.9894
## F-statistic: 2054 on 1 and 21 DF, p-value: < 2.2e-16length(wn2506_FLTL_lm$fitted.values)## [1] 23Scatterplot the residuals against the fitted values to check homogeneity of variance
plot(resid(wn2506_FLTL_lm) ~ fitted(wn2506_FLTL_lm))
Plot a histogram to check normality of residuals
qqnorm(wn2506_FLTL_lm$residuals)
qqline(wn2506_FLTL_lm$residuals)
fine
Fork length (FL_CM) - Standard length (SL_CM) relationship
wn2506_FLSL_plot <- ggplot(wndf_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 = "Pristipomoides aquilonaris FL-SL relationship (n = 24)",
fill = "Sex") +
scale_fill_viridis(option = "C",
discrete = TRUE,
begin = 0.2,
end = 0.8,
labels = c("N/A", "Female", "Male"))
wn2506_FLSL_plot
Linear model
wn2506_FLSL_lm <- lm(FL_CM ~ SL_CM, data = wndf_2506_LLrelationships)
print(summary(wn2506_FLSL_lm), digits = 4)##
## Call:
## lm(formula = FL_CM ~ SL_CM, data = wndf_2506_LLrelationships)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.67716 -0.16068 -0.00051 0.14302 0.74358
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.07725 0.50136 0.154 0.879
## SL_CM 1.10370 0.01765 62.540 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3588 on 22 degrees of freedom
## Multiple R-squared: 0.9944, Adjusted R-squared: 0.9942
## F-statistic: 3911 on 1 and 22 DF, p-value: < 2.2e-16length(wn2506_FLSL_lm$fitted.values)## [1] 24Scatterplot the residuals against the fitted values to check homogeneity of variance
plot(resid(wn2506_FLSL_lm) ~ fitted(wn2506_FLSL_lm))
Plot a histogram to check normality of residuals
qqnorm(resid(wn2506_FLSL_lm))
qqline(resid(wn2506_FLSL_lm))
Total length (TL_CM) - Standard length (SL_CM) relationship
wn2506_TLSL_plot <- ggplot(wndf_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 = "Pristipomoides aquilonaris TL-SL relationship (n = 23)",
fill = "Sex") +
scale_fill_viridis(option = "C",
discrete = TRUE,
begin = 0.2,
end = 0.8,
labels = c("N/A", "Female", "Male"))
wn2506_TLSL_plot
Linear model
wn2506_TLSL_lm <- lm(TL_CM ~ SL_CM, data = wndf_2506_LLrelationships)
summary(wn2506_TLSL_lm)##
## Call:
## lm(formula = TL_CM ~ SL_CM, data = wndf_2506_LLrelationships)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.81283 -0.30782 -0.03283 0.37343 1.23718
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.07279 0.79922 2.594 0.017 *
## SL_CM 1.17500 0.02781 42.244 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.51 on 21 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.9884, Adjusted R-squared: 0.9878
## F-statistic: 1785 on 1 and 21 DF, p-value: < 2.2e-16length(wn2506_TLSL_lm$fitted.values)## [1] 23Scatterplot the residuals against the fitted values to check homogeneity of variance
plot(resid(wn2506_TLSL_lm) ~ fitted(wn2506_TLSL_lm))
Plot a histogram to check normality of residuals
qqnorm(resid(wn2506_TLSL_lm))
qqline(resid(wn2506_TLSL_lm))
Fork length (FL_CM) - Weight (W_G) relationship
Load data file
wndf_2506_FL_W <- read.csv("wrn2506df.csv")
wndf_2506_FL_W$catch_month <- factor(wndf_2506_FL_W$catch_month)
wndf_2506_FL_W$year <- factor(wndf_2506_FL_W$year)
wndf_2506_FL_W$species <- factor(wndf_2506_FL_W$species)
wndf_2506_FL_W$common_name <- factor(wndf_2506_FL_W$common_name)
wndf_2506_FL_W$frozen <- factor(wndf_2506_FL_W$frozen)
wndf_2506_FL_W$catch_method <- factor(wndf_2506_FL_W$catch_method)
wndf_2506_FL_W$bait <- factor(wndf_2506_FL_W$bait)
wndf_2506_FL_W$area <- factor(wndf_2506_FL_W$area)
wndf_2506_FL_W$gutted <- factor(wndf_2506_FL_W$gutted)
wndf_2506_FL_W$sex <- factor(wndf_2506_FL_W$sex)
wndf_2506_FL_W$maturity <- factor(wndf_2506_FL_W$maturity)
# Remove gutted samples
wndf_2506_FL_W <- subset(filter(wndf_2506_FL_W, gutted == "N"))
# remove NA rows
wndf_2506_FL_W <- wndf_2506_FL_W%>%
drop_na(FL_CM)
View(wndf_2506_FL_W)
summary(wndf_2506_FL_W)## X catch_date catch_month landing_time
## Min. : 1.00 Length:24 Apr: 1 Length:24
## 1st Qu.: 6.75 Class :character Jun: 6 Class :character
## Median :12.50 Mode :character May: 1 Mode :character
## Mean :12.50 Sep:16
## 3rd Qu.:18.25
## Max. :24.00
##
## dissection_date dissection_time year month
## Length:24 Length:24 2023:16 Length:24
## Class :character Class :character 2024: 8 Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
## observer species common_name
## Length:24 Pristipomoides aquilonaris:24 Wenchman Snapper:24
## Class :character
## Mode :character
##
##
##
##
## fish_code frozen catch_method bait area
## Length:24 N: 8 Dropline :16 Japanese bait : 8 : 1
## Class :character Y:16 Redfish trap: 8 Japanese bait; Squid:16 B4:16
## Mode :character D4: 1
## D5: 6
##
##
##
## SL_CM FL_CM TL_CM HH_CM ED_CM
## Min. :19.50 Min. :21.60 Min. :25.50 Min. :4.30 Min. :1.500
## 1st Qu.:24.75 1st Qu.:27.27 1st Qu.:32.35 1st Qu.:5.05 1st Qu.:1.700
## Median :28.00 Median :31.15 Median :35.20 Median :5.70 Median :2.000
## Mean :28.10 Mean :31.10 Mean :35.53 Mean :5.60 Mean :2.057
## 3rd Qu.:31.43 3rd Qu.:34.60 3rd Qu.:38.70 3rd Qu.:6.05 3rd Qu.:2.450
## Max. :36.40 Max. :40.30 Max. :45.40 Max. :7.00 Max. :2.600
## NA's :1 NA's :17 NA's :17
## W_G finclip otoliths gutted
## Min. : 197.3 Length:24 Length:24 N:24
## 1st Qu.: 389.6 Class :character Class :character
## Median : 591.7 Mode :character Mode :character
## Mean : 626.0
## 3rd Qu.: 798.0
## Max. :1302.6
##
## gonads_present sex maturity G1L_CM G2L_CM
## Length:24 : 1 :16 Min. :1.500 Min. :2.100
## Class :character F:15 A: 3 1st Qu.:3.175 1st Qu.:3.350
## Mode :character M: 8 B: 5 Median :4.100 Median :3.900
## Mean :3.962 Mean :4.239
## 3rd Qu.:4.625 3rd Qu.:5.000
## Max. :7.000 Max. :7.100
## NA's :1
## GW_G stomach._sample
## Min. : 0.000 Length:24
## 1st Qu.: 0.550 Class :character
## Median : 0.850 Mode :character
## Mean : 3.008
## 3rd Qu.: 2.200
## Max. :35.300
## First look at the data with a scatterplot
fill = sex
ggplot(wndf_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 = "Pristipomoides aquilonaris, n = 24") +
scale_fill_viridis(option = "D", discrete = TRUE)
LM Linear model FL-W
wn2506_FLW_lm <- lm(FL_CM ~ W_G, data = wndf_2506_FL_W)
summary(wn2506_FLW_lm)##
## Call:
## lm(formula = FL_CM ~ W_G, data = wndf_2506_FL_W)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.4680 -0.5137 0.1927 0.7376 2.2354
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20.833723 0.605925 34.38 < 2e-16 ***
## W_G 0.016393 0.000888 18.46 7.05e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.181 on 22 degrees of freedom
## Multiple R-squared: 0.9394, Adjusted R-squared: 0.9366
## F-statistic: 340.8 on 1 and 22 DF, p-value: 7.05e-15Scatterplot the residuals against the fitted values to check homogeneity of variance
plot(resid(wn2506_FLW_lm) ~ fitted(wn2506_FLW_lm))
The scatter plot looks like it indicates a nonlinear relationship
Plot a histogram to check normality of residuals
hist(resid(wn2506_FLW_lm))
Plot the model
wn2506_FLW_lm_plottemp <- ggplot(wndf_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 = "Pristipomoides aquilonaris, n = 24") +
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
wn2506_FLW_lm_plot <- wn2506_FLW_lm_plottemp +
geom_text(x = 21, y = 1250,
label = expression(W~'(g)' == 0.016 * FL~'(cm)' + 20.834~"," ~ r^{2} == 0.937 ~ "(Linear)"),
colour = safe_colorblind_palette[1],
hjust = 0)
wn2506_FLW_lm_plot
NLS Models
NLS 1
wn2506_FLW_nls <- nls_multstart(W_G~A*FL_CM^B,
data = wndf_2506_FL_W,
iter = 500,
start_lower = c(A=0, B=2),
start_upper = c(A=0.5, B=3.5))
summary(wn2506_FLW_nls)##
## Formula: W_G ~ A * FL_CM^B
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## A 0.013261 0.005044 2.629 0.0153 *
## B 3.111114 0.107507 28.939 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 40.86 on 22 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
#check residuals
hist(resid(wn2506_FLW_nls))
Calculate R-squared
rss_fitted <- sum(residuals(wn2506_FLW_nls)^2)
rss_constant <- sum((wndf_2506_FL_W$W_G - mean(wndf_2506_FL_W$W_G))^2)
approx_r_squared <- 1 - rss_fitted / rss_constant
approx_r_squared## [1] 0.9792545Calculate AIC
AIC(wn2506_FLW_nls)## [1] 250.105Calculate the predicted W_kg values based on this NLS model
wn2506_FLW_nls_x <- c(21:41) #rounded actual min and max FL_CM values
wn2506_FLW_nls_predic <- data.frame(wn2506_FLW_nls_x)
wn2506_FLW_nls_f <- function(wn2506_FLW_nls_x) {
#formula using calculated parameters
return(0.013261*wn2506_FLW_nls_x^3.111114)
}
wn2506_FLW_nls_predic$wn2506_FLW_nls_y <- sapply(wn2506_FLW_nls_predic$wn2506_FLW_nls_x, wn2506_FLW_nls_f)
View(wn2506_FLW_nls_predic)Plot the model on the same graph as the other model
wn2506_FLW_nls_plottemp <- wn2506_FLW_lm_plot +
geom_smooth(aes(x = wn2506_FLW_nls_x, y = wn2506_FLW_nls_y),
data = wn2506_FLW_nls_predic,
method = "loess",
se = F,
color = safe_colorblind_palette[2],
linewidth = 0.75,
alpha = 0.75)
# add label with equation
wn2506_FLW_nls_plot <- wn2506_FLW_nls_plottemp +
geom_text(x = 21, y = 1150,
label = expression(W~'(g)' == 1.326%*%10^{-02} * FL~'(cm)'^{3.111}~"," ~ R^{2} == 0.979 ~ '(NLS)'),
colour = safe_colorblind_palette[2],
hjust = 0)
wn2506_FLW_nls_plot
Brms1: with nls parameters and a wide normal distribution around them as priors
set priors and run model
#set priors
wn2506_FLW_prior1 <- prior(normal(0.013261, 0.005) ,nlpar = "A")+prior(normal(3.111114, 1),nlpar = "B")
wn2506_FLW_brm1 <- brm(bf(W_G~A*FL_CM^B, A + B ~ 1, nl = TRUE),data = wndf_2506_FL_W,
prior = wn2506_FLW_prior1,
warmup = 1000,
iter = 5000,
chains = 4,
cores = 4,
control = list(adapt_delta = 0.8))
print(summary(wn2506_FLW_brm1), digits = 6)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: W_G ~ A * FL_CM^B
## A ~ 1
## B ~ 1
## Data: wndf_2506_FL_W (Number of observations: 24)
## 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.013994 0.003556 0.007703 0.021697 1.001664 2816 732
## B_Intercept 3.105263 0.074347 2.972295 3.265037 1.001890 2812 735
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 42.587197 6.771048 31.764878 57.633929 1.001334 3423 1037
##
## 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(wn2506_FLW_brm1))
Calculate R-squared
rss_fitted <- sum(residuals(wn2506_FLW_brm1)^2)
rss_constant <- sum((wndf_2506_FL_W$W_G - mean(wndf_2506_FL_W$W_G))^2)
approx_r_squared <- 1 - rss_fitted / rss_constant
approx_r_squared ## [1] 0.7013086smaller than the nls r-squared
Calculate the predicted W_G values based on this NLS model
wn2506_FLW_brm1_x <- c(21:41) #rounded actual min and max FL values
wn2506_FLW_brm1_predic <- data.frame(wn2506_FLW_brm1_x)
wn2506_FLW_brm1_f <- function(wn2506_FLW_brm1_x) {
#formula using calculated parameters
return(0.013720*wn2506_FLW_brm1_x^3.111616)
}
wn2506_FLW_brm1_predic$wn2506_FLW_brm1_y <- sapply(wn2506_FLW_brm1_predic$wn2506_FLW_brm1_x, wn2506_FLW_brm1_f)
View(wn2506_FLW_brm1_predic)Plot
wn2506_FLW_brm1_plottemp <- wn2506_FLW_nls_plot +
geom_smooth(aes(x = wn2506_FLW_brm1_x, y = wn2506_FLW_brm1_y),
data = wn2506_FLW_brm1_predic,
method = "loess",
se = F,
color = safe_colorblind_palette[4],
linewidth = 0.75,
alpha = 0.75)
# add label with equation
wn2506_FLW_brm1_plot <- wn2506_FLW_brm1_plottemp +
geom_text(x = 21, y = 1050,
label = expression(W~'(g)' == 1.372%*%10^{-02} * FL~'(cm)'^{3.112}~"," ~ R^{2} == 0.698 ~ '(BRMS)'),
colour = safe_colorblind_palette[4],
hjust = 0)
wn2506_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.021336-0.007416 ## [1] 0.01392#B UCI-LCI
3.275222-2.976599## [1] 0.298623set priors and run model
#set priors
wn2506_FLW_prior2 <- prior(normal(0.013720, 0.01392) ,nlpar = "A")+prior(normal(3.111616, 0.298623),nlpar = "B")
wn2506_FLW_brm2 <- brm(bf(W_G~A*FL_CM^B, A + B ~ 1, nl = TRUE),data = wndf_2506_FL_W,
prior = wn2506_FLW_prior2,
warmup = 1000,
iter = 5000,
chains = 4,
cores = 4,
control = list(adapt_delta = 0.8))
print(summary(wn2506_FLW_brm2), digits = 6)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: W_G ~ A * FL_CM^B
## A ~ 1
## B ~ 1
## Data: wndf_2506_FL_W (Number of observations: 24)
## 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.015504 0.005125 0.007755 0.027554 1.001621 3369 3477
## B_Intercept 3.081659 0.091811 2.903997 3.262570 1.001830 3376 3454
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 42.977106 6.807984 32.190457 58.549975 1.001465 5023 5490
##
## 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(wn2506_FLW_brm2))
Calculate R-squared
rss_fitted <- sum(residuals(wn2506_FLW_brm2)^2)
rss_constant <- sum((wndf_2506_FL_W$W_G - mean(wndf_2506_FL_W$W_G))^2)
approx_r_squared <- 1 - rss_fitted / rss_constant
approx_r_squared ## [1] 0.6964465Calculate the predicted W_G values based on this NLS model
wn2506_FLW_brm2_x <- c(21:41) #rounded actual min and max FL values
wn2506_FLW_brm2_predic <- data.frame(wn2506_FLW_brm2_x)
wn2506_FLW_brm2_f <- function(wn2506_FLW_brm2_x) {
#formula using calculated parameters
return(0.015489*wn2506_FLW_brm2_x^3.081884)
}
wn2506_FLW_brm2_predic$wn2506_FLW_brm2_y <- sapply(wn2506_FLW_brm2_predic$wn2506_FLW_brm2_x, wn2506_FLW_brm2_f)
View(wn2506_FLW_brm2_predic)Plot
wn2506_FLW_brm2_plottemp <- wn2506_FLW_brm1_plot +
geom_smooth(aes(x = wn2506_FLW_brm2_x, y = wn2506_FLW_brm2_y),
data = wn2506_FLW_brm2_predic,
method = "loess",
se = F,
color = safe_colorblind_palette[5],
linewidth = 0.75)
# add label with equation
wn2506_FLW_brm2_plot <- wn2506_FLW_brm2_plottemp +
geom_text(x = 21, y = 950,
label = expression(W~'(g)' == 1.549%*%10^{-02} * FL~'(cm)'^{3.082}~"," ~ R^{2} == 0.694 ~ '(BRMS)'),
colour = safe_colorblind_palette[5],
hjust = 0)
wn2506_FLW_brm2_plot
Interesting that the second brms is a worse fit than the first
##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
wn2506_FLW_nls_plot_title <- expression(italic("P. aquilonaris")~"FL-W Relationship (n = 24)")
wn2506_FLW_nls_plottemp <- ggplot(wndf_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 = wn2506_FLW_nls_plot_title) +
theme_classic() +
geom_smooth(aes(x = wn2506_FLW_nls_x, y = wn2506_FLW_nls_y),
data = wn2506_FLW_nls_predic,
method = "loess",
se = F,
color = turbo_colors_onboard[8],
linewidth = 0.75,
alpha = 0.75) +
ylim(0, 1600) +
xlim(12, 48)
# add label with equation
wn2506_FLW_nls_plot <- wn2506_FLW_nls_plottemp +
geom_label(x = 12, y = 1500,
label = expression(W~'(g)' == 1.326%*%10^{-02} * FL~'(cm)'^{3.111}~"," ~ R^{2} == 0.979 ~ '(NLS)'),
fill = "#ffb18e",
hjust = 0)
wn2506_FLW_nls_plot
Save plot
ggsave("wn2506_FLW_nls_plot.png", plot = wn2506_FLW_nls_plot, width = 7.2, height = 4, dpi = 400)## `geom_smooth()` using formula = 'y ~ x'## 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
wndf_2506_TL_W <- read.csv("wrn2506df.csv")
wndf_2506_TL_W$catch_month <- factor(wndf_2506_TL_W$catch_month)
wndf_2506_TL_W$year <- factor(wndf_2506_TL_W$year)
wndf_2506_TL_W$species <- factor(wndf_2506_TL_W$species)
wndf_2506_TL_W$common_name <- factor(wndf_2506_TL_W$common_name)
wndf_2506_TL_W$frozen <- factor(wndf_2506_TL_W$frozen)
wndf_2506_TL_W$catch_method <- factor(wndf_2506_TL_W$catch_method)
wndf_2506_TL_W$bait <- factor(wndf_2506_TL_W$bait)
wndf_2506_TL_W$area <- factor(wndf_2506_TL_W$area)
wndf_2506_TL_W$gutted <- factor(wndf_2506_TL_W$gutted)
wndf_2506_TL_W$sex <- factor(wndf_2506_TL_W$sex)
wndf_2506_TL_W$maturity <- factor(wndf_2506_TL_W$maturity)
# Remove gutted samples
wndf_2506_TL_W <- subset(filter(wndf_2506_TL_W, gutted == "N"))
# remove NA rows
wndf_2506_TL_W <- wndf_2506_TL_W%>%
drop_na(TL_CM)
View(wndf_2506_TL_W)
summary(wndf_2506_TL_W)## X catch_date catch_month landing_time
## Min. : 1.0 Length:23 Apr: 1 Length:23
## 1st Qu.: 6.5 Class :character Jun: 6 Class :character
## Median :12.0 Mode :character May: 1 Mode :character
## Mean :12.3 Sep:15
## 3rd Qu.:18.5
## Max. :24.0
##
## dissection_date dissection_time year month
## Length:23 Length:23 2023:15 Length:23
## Class :character Class :character 2024: 8 Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
## observer species common_name
## Length:23 Pristipomoides aquilonaris:23 Wenchman Snapper:23
## Class :character
## Mode :character
##
##
##
##
## fish_code frozen catch_method bait area
## Length:23 N: 8 Dropline :15 Japanese bait : 8 : 1
## Class :character Y:15 Redfish trap: 8 Japanese bait; Squid:15 B4:15
## Mode :character D4: 1
## D5: 6
##
##
##
## SL_CM FL_CM TL_CM HH_CM ED_CM
## Min. :20.20 Min. :22.10 Min. :25.50 Min. :4.30 Min. :1.500
## 1st Qu.:25.55 1st Qu.:28.25 1st Qu.:32.35 1st Qu.:5.05 1st Qu.:1.700
## Median :28.20 Median :31.20 Median :35.20 Median :5.70 Median :2.000
## Mean :28.48 Mean :31.51 Mean :35.53 Mean :5.60 Mean :2.057
## 3rd Qu.:31.45 3rd Qu.:34.60 3rd Qu.:38.70 3rd Qu.:6.05 3rd Qu.:2.450
## Max. :36.40 Max. :40.30 Max. :45.40 Max. :7.00 Max. :2.600
## NA's :16 NA's :16
## W_G finclip otoliths gutted
## Min. : 217.5 Length:23 Length:23 N:23
## 1st Qu.: 426.6 Class :character Class :character
## Median : 593.5 Mode :character Mode :character
## Mean : 644.6
## 3rd Qu.: 801.0
## Max. :1302.6
##
## gonads_present sex maturity G1L_CM G2L_CM
## Length:23 : 1 :15 Min. :1.500 Min. :2.100
## Class :character F:15 A: 3 1st Qu.:3.200 1st Qu.:3.450
## Mode :character M: 7 B: 5 Median :4.200 Median :3.900
## Mean :4.009 Mean :4.286
## 3rd Qu.:4.650 3rd Qu.:5.050
## Max. :7.000 Max. :7.100
## NA's :1
## GW_G stomach._sample
## Min. : 0.000 Length:23
## 1st Qu.: 0.600 Class :character
## Median : 0.900 Mode :character
## Mean : 3.139
## 3rd Qu.: 2.400
## Max. :35.300
## First look at the data with a scatterplot
fill = sex
ggplot(wndf_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 = "Pristipomoides aquilonaris, n = 23") +
scale_fill_viridis(option = "D", discrete = TRUE)
Linear model
wn2506_TLW_lm1 <- lm(TL_CM ~ W_G, data = wndf_2506_TL_W)
summary(wn2506_TLW_lm1)##
## Call:
## lm(formula = TL_CM ~ W_G, data = wndf_2506_TL_W)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.8418 -0.4811 0.1872 0.6371 1.9711
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.468e+01 5.663e-01 43.58 < 2e-16 ***
## W_G 1.684e-02 8.138e-04 20.69 1.9e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.022 on 21 degrees of freedom
## Multiple R-squared: 0.9532, Adjusted R-squared: 0.951
## F-statistic: 428.2 on 1 and 21 DF, p-value: 1.897e-15Scatterplot the residuals against the fitted values to check homogeneity of variance
plot(resid(wn2506_TLW_lm1) ~ fitted(wn2506_TLW_lm1))
looks (slightly) nonlinear
Plot a histogram to check normality of residuals
qqnorm(resid(wn2506_TLW_lm1))
qqline(resid(wn2506_TLW_lm1))
Plot the model
wn2506_TLW_lm1_plottemp <- ggplot(wndf_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 = 23") +
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
wn2506_TLW_lm1_plot <- wn2506_TLW_lm1_plottemp +
geom_text(x = 25, y = 1300,
label = expression(W~'(g)' == 0.022 * TL~'(cm)' + 20.84~"," ~ r^{2} == 0.9118 ~ "(Linear)"),
colour = safe_colorblind_palette[1],
hjust = 0)
wn2506_TLW_lm1_plot
NLS
wn2506_TLW_nls <- nls_multstart(W_G~A*TL_CM^B,
data = wndf_2506_TL_W,
iter = 500,
start_lower = c(A=0, B=2),
start_upper = c(A=0.5, B=3.5))
summary(wn2506_TLW_nls)##
## Formula: W_G ~ A * TL_CM^B
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## A 0.006048 0.002236 2.705 0.0133 *
## B 3.226819 0.101106 31.915 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 36.98 on 21 degrees of freedom
##
## Number of iterations to convergence: 8
## Achieved convergence tolerance: 1.49e-08Plot a histogram of the residuals to check normality of residuals
qqnorm(resid(wn2506_TLW_nls))
qqline(resid(wn2506_TLW_nls))
Calculate R-squared
rss_fitted <- sum(residuals(wn2506_TLW_nls)^2)
rss_constant <- sum((wndf_2506_TL_W$W_G - mean(wndf_2506_TL_W$W_G))^2)
approx_r_squared <- 1 - rss_fitted / rss_constant
approx_r_squared## [1] 0.9818043Calculate AIC
AIC(wn2506_TLW_nls)## [1] 235.2593Calculate the predicted W_kg values based on this NLS model
wn2506_TLW_nls_x <- c(25:46) #rounded actual min and max TL_CM values
wn2506_TLW_nls_predic <- data.frame(wn2506_TLW_nls_x)
wn2506_TLW_nls_f <- function(wn2506_TLW_nls_x) {
#formula using calculated parameters
return(0.006048*wn2506_TLW_nls_x^3.226819)
}
wn2506_TLW_nls_predic$wn2506_TLW_nls_y <- sapply(wn2506_TLW_nls_predic$wn2506_TLW_nls_x, wn2506_TLW_nls_f)
View(wn2506_TLW_nls_predic)Plot the model on the same graph as the other model
wn2506_TLW_nls_plottemp <- wn2506_TLW_lm1_plot +
geom_smooth(aes(x = wn2506_TLW_nls_x, y = wn2506_TLW_nls_y),
data = wn2506_TLW_nls_predic,
method = "loess",
se = F,
color = safe_colorblind_palette[2],
linewidth = 0.75)
# add label with equation
wn2506_TLW_nls_plot <- wn2506_TLW_nls_plottemp +
geom_text(x = 25, y = 1200,
label = expression(W~'(g)' == 6.048%*%10^{-03} * TL~'(cm)'^{3.227}~"," ~ R^{2} == 0.982 ~ '(NLS)'),
colour = safe_colorblind_palette[2],
hjust = 0)
wn2506_TLW_nls_plot
Brms: with nls parameters to start
set priors and run model
#set priors
wn2506_TLW_prior1 <- prior(normal(0.006048, 0.002236) ,nlpar = "A")+prior(normal(3.226819, 0.101106),nlpar = "B")
wn2506_TLW_brm1 <- brm(bf(W_G~A*TL_CM^B, A + B ~ 1, nl = TRUE),data = wndf_2506_TL_W,
prior = wn2506_TLW_prior1,
warmup = 1000,
iter = 5000,
chains = 4,
cores = 4,
control = list(adapt_delta = 0.8))
print(summary(wn2506_TLW_brm1), digits = 6)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: W_G ~ A * TL_CM^B
## A ~ 1
## B ~ 1
## Data: wndf_2506_TL_W (Number of observations: 23)
## 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.006284 0.001305 0.003959 0.009036 1.000777 4246 3380
## B_Intercept 3.222253 0.057649 3.117226 3.342151 1.000776 4229 3380
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 38.639286 6.198718 28.863958 53.168977 0.999957 5812 5060
##
## 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(wn2506_TLW_brm1))
Calculate R-squared
rss_fitted <- sum(residuals(wn2506_TLW_brm1)^2)
rss_constant <- sum((wndf_2506_TL_W$W_G - mean(wndf_2506_TL_W$W_G))^2)
approx_r_squared <- 1 - rss_fitted / rss_constant
approx_r_squared ## [1] 0.7384464Calculate the predicted W_G values based on this NLS model
wn2506_TLW_brm1_x <- c(25:46) #rounded actual min and max TL values
wn2506_TLW_brm1_predic <- data.frame(wn2506_TLW_brm1_x)
wn2506_TLW_brm1_f <- function(wn2506_TLW_brm1_x) {
#formula using calculated parameters
return(0.006336*wn2506_TLW_brm1_x^3.219897)
}
wn2506_TLW_brm1_predic$wn2506_TLW_brm1_y <- sapply(wn2506_TLW_brm1_predic$wn2506_TLW_brm1_x, wn2506_TLW_brm1_f)
View(wn2506_TLW_brm1_predic)Plot
wn2506_TLW_brm1_plottemp <- wn2506_TLW_nls_plot +
geom_smooth(aes(x = wn2506_TLW_brm1_x, y = wn2506_TLW_brm1_y),
data = wn2506_TLW_brm1_predic,
method = "loess",
se = F,
color = safe_colorblind_palette[4],
linewidth = 0.75)
# add label with equation
wn2506_TLW_brm1_plot <- wn2506_TLW_brm1_plottemp +
geom_text(x = 25, y = 1100,
label = expression(W~'(g)' == 6.336%*%10^{-03} * TL~'(cm)'^{"3.220"}~"," ~ R^{2} == 0.734 ~ '(BRMS)'),
colour = safe_colorblind_palette[4],
hjust = 0)
wn2506_TLW_brm1_plot
The NLS model and the BRMS models are very similar which is a good sign
##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
wndf_2506_SL_W <- read.csv("wrn2506df.csv")
wndf_2506_SL_W$catch_month <- factor(wndf_2506_SL_W$catch_month)
wndf_2506_SL_W$year <- factor(wndf_2506_SL_W$year)
wndf_2506_SL_W$species <- factor(wndf_2506_SL_W$species)
wndf_2506_SL_W$common_name <- factor(wndf_2506_SL_W$common_name)
wndf_2506_SL_W$frozen <- factor(wndf_2506_SL_W$frozen)
wndf_2506_SL_W$catch_method <- factor(wndf_2506_SL_W$catch_method)
wndf_2506_SL_W$bait <- factor(wndf_2506_SL_W$bait)
wndf_2506_SL_W$area <- factor(wndf_2506_SL_W$area)
wndf_2506_SL_W$gutted <- factor(wndf_2506_SL_W$gutted)
wndf_2506_SL_W$sex <- factor(wndf_2506_SL_W$sex)
wndf_2506_SL_W$maturity <- factor(wndf_2506_SL_W$maturity)
# Remove gutted samples
wndf_2506_SL_W <- subset(filter(wndf_2506_SL_W, gutted == "N"))
# remove NA rows
wndf_2506_SL_W <- wndf_2506_SL_W%>%
drop_na(SL_CM)
View(wndf_2506_SL_W)
summary(wndf_2506_SL_W)## X catch_date catch_month landing_time
## Min. : 1.00 Length:24 Apr: 1 Length:24
## 1st Qu.: 6.75 Class :character Jun: 6 Class :character
## Median :12.50 Mode :character May: 1 Mode :character
## Mean :12.50 Sep:16
## 3rd Qu.:18.25
## Max. :24.00
##
## dissection_date dissection_time year month
## Length:24 Length:24 2023:16 Length:24
## Class :character Class :character 2024: 8 Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
## observer species common_name
## Length:24 Pristipomoides aquilonaris:24 Wenchman Snapper:24
## Class :character
## Mode :character
##
##
##
##
## fish_code frozen catch_method bait area
## Length:24 N: 8 Dropline :16 Japanese bait : 8 : 1
## Class :character Y:16 Redfish trap: 8 Japanese bait; Squid:16 B4:16
## Mode :character D4: 1
## D5: 6
##
##
##
## SL_CM FL_CM TL_CM HH_CM ED_CM
## Min. :19.50 Min. :21.60 Min. :25.50 Min. :4.30 Min. :1.500
## 1st Qu.:24.75 1st Qu.:27.27 1st Qu.:32.35 1st Qu.:5.05 1st Qu.:1.700
## Median :28.00 Median :31.15 Median :35.20 Median :5.70 Median :2.000
## Mean :28.10 Mean :31.10 Mean :35.53 Mean :5.60 Mean :2.057
## 3rd Qu.:31.43 3rd Qu.:34.60 3rd Qu.:38.70 3rd Qu.:6.05 3rd Qu.:2.450
## Max. :36.40 Max. :40.30 Max. :45.40 Max. :7.00 Max. :2.600
## NA's :1 NA's :17 NA's :17
## W_G finclip otoliths gutted
## Min. : 197.3 Length:24 Length:24 N:24
## 1st Qu.: 389.6 Class :character Class :character
## Median : 591.7 Mode :character Mode :character
## Mean : 626.0
## 3rd Qu.: 798.0
## Max. :1302.6
##
## gonads_present sex maturity G1L_CM G2L_CM
## Length:24 : 1 :16 Min. :1.500 Min. :2.100
## Class :character F:15 A: 3 1st Qu.:3.175 1st Qu.:3.350
## Mode :character M: 8 B: 5 Median :4.100 Median :3.900
## Mean :3.962 Mean :4.239
## 3rd Qu.:4.625 3rd Qu.:5.000
## Max. :7.000 Max. :7.100
## NA's :1
## GW_G stomach._sample
## Min. : 0.000 Length:24
## 1st Qu.: 0.550 Class :character
## Median : 0.850 Mode :character
## Mean : 3.008
## 3rd Qu.: 2.200
## Max. :35.300
## First look at the data with a scatterplot
fill = sex
ggplot(wndf_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 = "Pristipomoides aquilonaris, n = 24") +
scale_fill_viridis(option = "D", discrete = TRUE)
LM Linear model FL-W
wn2506_FLW_lm <- lm(SL_CM ~ W_G, data = wndf_2506_SL_W)
summary(wn2506_FLW_lm)##
## Call:
## lm(formula = SL_CM ~ W_G, data = wndf_2506_SL_W)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.22618 -0.31296 0.06754 0.67684 1.42075
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.879e+01 5.081e-01 36.98 < 2e-16 ***
## W_G 1.488e-02 7.446e-04 19.98 1.36e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9907 on 22 degrees of freedom
## Multiple R-squared: 0.9478, Adjusted R-squared: 0.9454
## F-statistic: 399.2 on 1 and 22 DF, p-value: 1.36e-15Scatterplot the residuals against the fitted values to check homogeneity of variance
plot(resid(wn2506_FLW_lm) ~ fitted(wn2506_FLW_lm))
looks like a nonlinear relationship
Plot a histogram to check normality of residuals
qqnorm(resid(wn2506_FLW_lm))
qqline(resid(wn2506_FLW_lm))
not normal
Plot the model
wn2506_SLW_lm_plottemp <- ggplot(wndf_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 = "Pristipomoides aquilonaris, n = 24") +
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
wn2506_SLW_lm_plot <- wn2506_SLW_lm_plottemp +
geom_text(x = 19, y = 1300,
label = expression(W~'(g)' == 0.029 * SL~'(cm)' + 13.611~"," ~ r^{2} == 0.9529 ~ "(Linear)"),
colour = safe_colorblind_palette[1],
hjust = 0)
wn2506_SLW_lm_plot
NLS Models
NLS 1
wn2506_SLW_nls <- nls_multstart(W_G~A*SL_CM^B,
data = wndf_2506_SL_W,
iter = 500,
start_lower = c(A=0, B=2),
start_upper = c(A=0.5, B=3.5))
print(summary(wn2506_SLW_nls), digits = 4)##
## Formula: W_G ~ A * SL_CM^B
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## A 0.020177 0.006796 2.969 0.00708 **
## B 3.080229 0.097978 31.438 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 37.81 on 22 degrees of freedom
##
## Number of iterations to convergence: 5
## Achieved convergence tolerance: 1.49e-08Plot a histogram of the residuals to check normality of residuals
#check residuals
hist(resid(wn2506_SLW_nls))
Calculate R-squared
rss_fitted <- sum(residuals(wn2506_SLW_nls)^2)
rss_constant <- sum((wndf_2506_SL_W$W_G - mean(wndf_2506_SL_W$W_G))^2)
approx_r_squared <- 1 - rss_fitted / rss_constant
approx_r_squared## [1] 0.9822326Calculate AIC
AIC(wn2506_SLW_nls)## [1] 246.3858Calculate the predicted W_g values based on this NLS model
wn2506_SLW_nls_x <- c(19:37) #rounded actual min and max FL_CM values
wn2506_SLW_nls_predic <- data.frame(wn2506_SLW_nls_x)
wn2506_SLW_nls_f <- function(wn2506_SLW_nls_x) {
#formula using calculated parameters
return(0.020177*wn2506_SLW_nls_x^3.080229)
}
wn2506_SLW_nls_predic$wn2506_SLW_nls_y <- sapply(wn2506_SLW_nls_predic$wn2506_SLW_nls_x, wn2506_SLW_nls_f)
View(wn2506_SLW_nls_predic)Plot the model on the same graph as the other model
wn2506_SLW_nls_plottemp <- wn2506_SLW_lm_plot +
geom_smooth(aes(x = wn2506_SLW_nls_x, y = wn2506_SLW_nls_y),
data = wn2506_SLW_nls_predic,
method = "loess",
se = F,
color = safe_colorblind_palette[2],
linewidth = 0.75,
alpha = 0.75)
# add label with equation
wn2506_SLW_nls_plot <- wn2506_SLW_nls_plottemp +
geom_text(x = 19, y = 1200,
label = expression(W~'(g)' == 2.018%*%10^{-02} * SL~'(cm)'^{"3.080"}~"," ~ R^{2} == 0.982 ~ '(NLS)'),
colour = safe_colorblind_palette[2],
hjust = 0)
wn2506_SLW_nls_plot
Brms: with nls parameters to start
set priors and run model
#set priors
wn2506_SLW_prior1 <- prior(normal(0.020177, 0.006796) ,nlpar = "A")+prior(normal(3.080229, 0.097978),nlpar = "B")
wn2506_SLW_brm1 <- brm(bf(W_G~A*SL_CM^B, A + B ~ 1, nl = TRUE),data = wndf_2506_SL_W,
prior = wn2506_SLW_prior1,
warmup = 1000,
iter = 5000,
chains = 4,
cores = 4,
control = list(adapt_delta = 0.8))
print(summary(wn2506_SLW_brm1), digits = 6)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: W_G ~ A * SL_CM^B
## A ~ 1
## B ~ 1
## Data: wndf_2506_SL_W (Number of observations: 24)
## 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.020924 0.003975 0.013937 0.029242 1.001445 4246 5245
## B_Intercept 3.074816 0.055568 2.971867 3.187167 1.001485 4219 5295
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 39.462358 6.307604 29.318102 54.180032 1.000211 5787 5393
##
## 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(wn2506_SLW_brm1))
Calculate R-squared
rss_fitted <- sum(residuals(wn2506_SLW_brm1)^2)
rss_constant <- sum((wndf_2506_SL_W$W_G - mean(wndf_2506_SL_W$W_G))^2)
approx_r_squared <- 1 - rss_fitted / rss_constant
approx_r_squared ## [1] 0.7469714Calculate the predicted W_G values based on this NLS model
wn2506_SLW_brm1_x <- c(19:37) #rounded actual min and max SL values
wn2506_SLW_brm1_predic <- data.frame(wn2506_SLW_brm1_x)
wn2506_SLW_brm1_f <- function(wn2506_SLW_brm1_x) {
#formula using calculated parameters
return(0.020888*wn2506_SLW_brm1_x^3.075490)
}
wn2506_SLW_brm1_predic$wn2506_SLW_brm1_y <- sapply(wn2506_SLW_brm1_predic$wn2506_SLW_brm1_x, wn2506_SLW_brm1_f)
View(wn2506_SLW_brm1_predic)Plot
wn2506_SLW_brm1_plottemp <- wn2506_SLW_nls_plot +
geom_smooth(aes(x = wn2506_SLW_brm1_x, y = wn2506_SLW_brm1_y),
data = wn2506_SLW_brm1_predic,
method = "loess",
se = F,
color = safe_colorblind_palette[3],
linewidth = 0.75)
# add label with equation
wn2506_SLW_brm1_plot <- wn2506_SLW_brm1_plottemp +
geom_text(x = 19, y = 1100,
label = expression(W~'(g)' == 2.068%*%10^{-02} * SL~'(cm)'^{3.077}~"," ~ R^{2} == 0.746 ~ '(BRMS)'),
colour = safe_colorblind_palette[3],
hjust = 0)
wn2506_SLW_brm1_plot
Final decision: NLS
Maturity
Load data file
wndf_2506_maturity <- read.csv("wrn2506df.csv")
wndf_2506_maturity$catch_month <- factor(wndf_2506_maturity$catch_month)
wndf_2506_maturity$year <- factor(wndf_2506_maturity$year)
wndf_2506_maturity$species <- factor(wndf_2506_maturity$species)
wndf_2506_maturity$common_name <- factor(wndf_2506_maturity$common_name)
wndf_2506_maturity$frozen <- factor(wndf_2506_maturity$frozen)
wndf_2506_maturity$catch_method <- factor(wndf_2506_maturity$catch_method)
wndf_2506_maturity$bait <- factor(wndf_2506_maturity$bait)
wndf_2506_maturity$area <- factor(wndf_2506_maturity$area)
wndf_2506_maturity$gutted <- factor(wndf_2506_maturity$gutted)
wndf_2506_maturity$sex <- factor(wndf_2506_maturity$sex)
wndf_2506_maturity$maturity <- factor(wndf_2506_maturity$maturity)
#drop the columns with neither sex nor maturity measurements
wndf_2506_maturity <- wndf_2506_maturity %>%
filter(maturity != "" | sex != "")calculate sex ratio
wn_totalfemale <- wndf_2506_maturity %>%
filter(sex == "F") %>%
count()
wn_totalmale <- wndf_2506_maturity %>%
filter(sex == "M") %>%
count()
wn_sexratio <- wn_totalfemale/wn_totalmale
wn_sexratio #female number## n
## 1 1.875Descriptive statistics
wn_sex_summary <- wndf_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)
)
wn_sex_summary## # A tibble: 3 × 8
## sex count mean sd median min max IQR
## <fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 "" 1 27 NA 27 27 27 0
## 2 "F" 15 31.9 4.91 33.1 22.1 40.3 6.8
## 3 "M" 8 30.1 4.38 30.8 21.6 35.2 3.62wndf_2506_maturity <- wndf_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))
wn_maturity_summary <- wndf_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)
)
wn_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 3 28.3 2.46 27 26.7 31.1 2.2
## 2 mature 5 27.7 4.23 27.2 22.1 33.1 4.6
## 3 <NA> 16 32.7 4.48 33.8 21.6 40.3 4.88