R. aurorubens body size models
Naomi Bakken
2024-08-15
This R markdown file includes all the R code used to analyze Vermilion Snapper (Rhomboplites aurorubens) 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
vsdf_2506_LLrelationships <- read.csv("vs2506df.csv")
vsdf_2506_LLrelationships$catch_month <- factor(vsdf_2506_LLrelationships$catch_month)
vsdf_2506_LLrelationships$year <- factor(vsdf_2506_LLrelationships$year)
vsdf_2506_LLrelationships$species <- factor(vsdf_2506_LLrelationships$species)
vsdf_2506_LLrelationships$common_name <- factor(vsdf_2506_LLrelationships$common_name)
vsdf_2506_LLrelationships$frozen <- factor(vsdf_2506_LLrelationships$frozen)
vsdf_2506_LLrelationships$catch_method <- factor(vsdf_2506_LLrelationships$catch_method)
vsdf_2506_LLrelationships$bait <- factor(vsdf_2506_LLrelationships$bait)
vsdf_2506_LLrelationships$area <- factor(vsdf_2506_LLrelationships$area)
vsdf_2506_LLrelationships$gutted <- factor(vsdf_2506_LLrelationships$gutted)
vsdf_2506_LLrelationships$sex <- factor(vsdf_2506_LLrelationships$sex)
vsdf_2506_LLrelationships$maturity <- factor(vsdf_2506_LLrelationships$maturity)
summary(vsdf_2506_LLrelationships)## X catch_date catch_month landing_time
## Min. : 1 Length:61 Apr: 1 Length:61
## 1st Qu.:16 Class :character Jul: 9 Class :character
## Median :31 Mode :character Jun:15 Mode :character
## Mean :31 May: 5
## 3rd Qu.:46 Sep:31
## Max. :61
##
## dissection_date dissection_time year month
## Length:61 Length:61 2023:49 Length:61
## Class :character Class :character 2024:12 Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
## observer species common_name
## Length:61 Rhomboplites aurorubens:61 Vermilion Snapper:61
## Class :character
## Mode :character
##
##
##
##
## fish_code frozen catch_method bait area
## Length:61 N:12 :49 :49 :19
## Class :character Y:49 Redfish trap:12 Japanese bait:12 C2:30
## Mode :character D3: 2
## D4: 1
## D5: 9
##
##
## SL_CM FL_CM TL_CM HH_CM
## Min. :16.90 Min. :14.00 Min. :14.90 Min. :3.400
## 1st Qu.:18.70 1st Qu.:16.50 1st Qu.:17.90 1st Qu.:3.700
## Median :19.40 Median :19.90 Median :21.60 Median :4.000
## Mean :20.98 Mean :20.76 Mean :22.81 Mean :4.245
## 3rd Qu.:24.00 3rd Qu.:25.20 3rd Qu.:28.00 3rd Qu.:4.800
## Max. :27.00 Max. :29.90 Max. :34.50 Max. :5.600
## NA's :48 NA's :50
## ED_CM W_G finclip otoliths gutted
## Min. :1.100 Min. : 52.7 Length:61 Length:61 N:61
## 1st Qu.:1.500 1st Qu.: 81.4 Class :character Class :character
## Median :1.600 Median :132.2 Mode :character Mode :character
## Mean :1.618 Mean :182.6
## 3rd Qu.:1.800 3rd Qu.:281.2
## Max. :2.100 Max. :481.6
## NA's :50
## gonads_present sex maturity G1L_CM G2L_CM
## Length:61 :12 :49 Min. :1.000 Min. :0.500
## Class :character F:33 B: 4 1st Qu.:2.500 1st Qu.:2.300
## Mode :character M:16 C: 8 Median :3.500 Median :3.200
## Mean :3.668 Mean :3.532
## 3rd Qu.:4.600 3rd Qu.:4.800
## Max. :6.500 Max. :6.500
## NA's :8 NA's :8
## GW_G stomach._sample
## Length:61 Length:61
## Class :character Class :character
## Mode :character Mode :character
##
##
##
## View(vsdf_2506_LLrelationships)
sumtable(vsdf_2506_LLrelationships, digits = 4)| Variable | N | Mean | Std. Dev. | Min | Pctl. 25 | Pctl. 75 | Max |
|---|---|---|---|---|---|---|---|
| X | 61 | 31 | 17.75 | 1 | 16 | 46 | 61 |
| catch_month | 61 | ||||||
| … Apr | 1 | 1.64% | |||||
| … Jul | 9 | 14.75% | |||||
| … Jun | 15 | 24.59% | |||||
| … May | 5 | 8.2% | |||||
| … Sep | 31 | 50.82% | |||||
| landing_time | 61 | ||||||
| … | 49 | 80.33% | |||||
| … 11:00 | 1 | 1.64% | |||||
| … 12:45 | 1 | 1.64% | |||||
| … 13:30 | 1 | 1.64% | |||||
| … 14:10 | 4 | 6.56% | |||||
| … 14:45 | 5 | 8.2% | |||||
| year | 61 | ||||||
| … 2023 | 49 | 80.33% | |||||
| … 2024 | 12 | 19.67% | |||||
| observer | 61 | ||||||
| … HK | 17 | 27.87% | |||||
| … LJ | 27 | 44.26% | |||||
| … NB | 5 | 8.2% | |||||
| … QH | 5 | 8.2% | |||||
| … SH | 1 | 1.64% | |||||
| … SM | 6 | 9.84% | |||||
| species | 61 | ||||||
| … Rhomboplites aurorubens | 61 | 100% | |||||
| common_name | 61 | ||||||
| … Vermilion Snapper | 61 | 100% | |||||
| frozen | 61 | ||||||
| … N | 12 | 19.67% | |||||
| … Y | 49 | 80.33% | |||||
| catch_method | 61 | ||||||
| … | 49 | 80.33% | |||||
| … Redfish trap | 12 | 19.67% | |||||
| bait | 61 | ||||||
| … | 49 | 80.33% | |||||
| … Japanese bait | 12 | 19.67% | |||||
| area | 61 | ||||||
| … | 19 | 31.15% | |||||
| … C2 | 30 | 49.18% | |||||
| … D3 | 2 | 3.28% | |||||
| … D4 | 1 | 1.64% | |||||
| … D5 | 9 | 14.75% | |||||
| SL_CM | 13 | 20.98 | 3.549 | 16.9 | 18.7 | 24 | 27 |
| FL_CM | 61 | 20.76 | 4.812 | 14 | 16.5 | 25.2 | 29.9 |
| TL_CM | 61 | 22.81 | 5.584 | 14.9 | 17.9 | 28 | 34.5 |
| HH_CM | 11 | 4.245 | 0.716 | 3.4 | 3.7 | 4.8 | 5.6 |
| ED_CM | 11 | 1.618 | 0.275 | 1.1 | 1.5 | 1.8 | 2.1 |
| W_G | 61 | 182.6 | 117.9 | 52.7 | 81.4 | 281.2 | 481.6 |
| finclip | 61 | ||||||
| … | 29 | 47.54% | |||||
| … N | 3 | 4.92% | |||||
| … Y | 29 | 47.54% | |||||
| otoliths | 61 | ||||||
| … Y | 61 | 100% | |||||
| gutted | 61 | ||||||
| … N | 61 | 100% | |||||
| gonads_present | 61 | ||||||
| … N | 5 | 8.2% | |||||
| … Y | 56 | 91.8% | |||||
| sex | 61 | ||||||
| … | 12 | 19.67% | |||||
| … F | 33 | 54.1% | |||||
| … M | 16 | 26.23% | |||||
| maturity | 61 | ||||||
| … | 49 | 80.33% | |||||
| … B | 4 | 6.56% | |||||
| … C | 8 | 13.11% | |||||
| G1L_CM | 53 | 3.668 | 1.327 | 1 | 2.5 | 4.6 | 6.5 |
| G2L_CM | 53 | 3.532 | 1.447 | 0.5 | 2.3 | 4.8 | 6.5 |
| stomach._sample | 61 | ||||||
| … N | 61 | 100% |
Fork length (FL_CM) - Total length (TL_CM) relationship
vs2506_FLTL_plotbase <- ggplot(vsdf_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 = "Rhomboplites aurorubens FL-TL relationship (n = 61)",
fill = "Sex") +
scale_fill_viridis(option = "C",
discrete = TRUE,
begin = 0.2,
end = 0.8,
labels = c("N/A", "Female", "Male"))
vs2506_FLTL_plotbase
Linear model
vs2506_FLTL_lm <- lm(FL_CM ~ TL_CM, data = vsdf_2506_LLrelationships)
summary(vs2506_FLTL_lm)##
## Call:
## lm(formula = FL_CM ~ TL_CM, data = vsdf_2506_LLrelationships)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.1749 -0.2163 0.0028 0.3304 0.8571
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.29507 0.36784 3.521 0.000837 ***
## TL_CM 0.85338 0.01567 54.452 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6779 on 59 degrees of freedom
## Multiple R-squared: 0.9805, Adjusted R-squared: 0.9802
## F-statistic: 2965 on 1 and 59 DF, p-value: < 2.2e-16length(vs2506_FLTL_lm$fitted.values)## [1] 61Scatterplot the residuals against the fitted values to check homogeneity of variance
plot(resid(vs2506_FLTL_lm) ~ fitted(vs2506_FLTL_lm))
fine other than one clear outlier
Plot a histogram to check normality of residuals
qqnorm(vs2506_FLTL_lm$residuals)
qqline(vs2506_FLTL_lm$residuals)
same
Fork length (FL_CM) - Standard length (SL_CM) relationship
vs2506_FLSL_plot <- ggplot(vsdf_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 = "Rhomboplites aurorubens FL-SL relationship (n = 13)",
fill = "Sex") +
scale_fill_viridis(option = "C",
discrete = TRUE,
begin = 0.2,
end = 0.8,
labels = c("N/A", "Female", "Male"))
vs2506_FLSL_plot
Linear model
vs2506_FLSL_lm <- lm(FL_CM ~ SL_CM, data = vsdf_2506_LLrelationships)
print(summary(vs2506_FLSL_lm), digits = 4)##
## Call:
## lm(formula = FL_CM ~ SL_CM, data = vsdf_2506_LLrelationships)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.7771 -0.3713 -0.1872 0.3876 0.8901
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.60205 0.98279 2.648 0.0227 *
## SL_CM 1.01933 0.04623 22.050 1.87e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5683 on 11 degrees of freedom
## (48 observations deleted due to missingness)
## Multiple R-squared: 0.9779, Adjusted R-squared: 0.9759
## F-statistic: 486.2 on 1 and 11 DF, p-value: 1.871e-10length(vs2506_FLSL_lm$fitted.values)## [1] 13Scatterplot the residuals against the fitted values to check homogeneity of variance
plot(resid(vs2506_FLSL_lm) ~ fitted(vs2506_FLSL_lm))
Plot a histogram to check normality of residuals
qqnorm(resid(vs2506_FLSL_lm))
qqline(resid(vs2506_FLSL_lm))
Total length (TL_CM) - Standard length (SL_CM) relationship
vs2506_TLSL_plot <- ggplot(vsdf_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 = "Rhomboplites aurorubens TL-SL relationship (n = 13)",
fill = "Sex") +
scale_fill_viridis(option = "C",
discrete = TRUE,
begin = 0.2,
end = 0.8,
labels = c("N/A", "Female", "Male"))
vs2506_TLSL_plot
Linear model
vs2506_TLSL_lm <- lm(TL_CM ~ SL_CM, data = vsdf_2506_LLrelationships)
print(summary(vs2506_TLSL_lm), digits = 5)##
## Call:
## lm(formula = TL_CM ~ SL_CM, data = vsdf_2506_LLrelationships)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.78081 -0.56176 -0.31703 0.09352 4.02281
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.60899 2.45810 1.0614 0.3113
## SL_CM 1.14803 0.11562 9.9292 7.937e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.4214 on 11 degrees of freedom
## (48 observations deleted due to missingness)
## Multiple R-squared: 0.89963, Adjusted R-squared: 0.8905
## F-statistic: 98.589 on 1 and 11 DF, p-value: 7.9365e-07length(vs2506_TLSL_lm$fitted.values)## [1] 13Scatterplot the residuals against the fitted values to check homogeneity of variance
plot(resid(vs2506_TLSL_lm) ~ fitted(vs2506_TLSL_lm))
Plot a histogram to check normality of residuals
qqnorm(resid(vs2506_TLSL_lm))
qqline(resid(vs2506_TLSL_lm))
Fork length (FL_CM) - Weight (W_G) relationship
Load data file
vsdf_2506_FL_W <- read.csv("vs2506df.csv")
vsdf_2506_FL_W$catch_month <- factor(vsdf_2506_FL_W$catch_month)
vsdf_2506_FL_W$year <- factor(vsdf_2506_FL_W$year)
vsdf_2506_FL_W$species <- factor(vsdf_2506_FL_W$species)
vsdf_2506_FL_W$common_name <- factor(vsdf_2506_FL_W$common_name)
vsdf_2506_FL_W$frozen <- factor(vsdf_2506_FL_W$frozen)
vsdf_2506_FL_W$catch_method <- factor(vsdf_2506_FL_W$catch_method)
vsdf_2506_FL_W$bait <- factor(vsdf_2506_FL_W$bait)
vsdf_2506_FL_W$area <- factor(vsdf_2506_FL_W$area)
vsdf_2506_FL_W$gutted <- factor(vsdf_2506_FL_W$gutted)
vsdf_2506_FL_W$sex <- factor(vsdf_2506_FL_W$sex)
vsdf_2506_FL_W$maturity <- factor(vsdf_2506_FL_W$maturity)
# Remove gutted samples
vsdf_2506_FL_W <- subset(filter(vsdf_2506_FL_W, gutted == "N"))
# remove NA rows
vsdf_2506_FL_W <- vsdf_2506_FL_W%>%
drop_na(FL_CM)
View(vsdf_2506_FL_W)
summary(vsdf_2506_FL_W)## X catch_date catch_month landing_time
## Min. : 1 Length:61 Apr: 1 Length:61
## 1st Qu.:16 Class :character Jul: 9 Class :character
## Median :31 Mode :character Jun:15 Mode :character
## Mean :31 May: 5
## 3rd Qu.:46 Sep:31
## Max. :61
##
## dissection_date dissection_time year month
## Length:61 Length:61 2023:49 Length:61
## Class :character Class :character 2024:12 Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
## observer species common_name
## Length:61 Rhomboplites aurorubens:61 Vermilion Snapper:61
## Class :character
## Mode :character
##
##
##
##
## fish_code frozen catch_method bait area
## Length:61 N:12 :49 :49 :19
## Class :character Y:49 Redfish trap:12 Japanese bait:12 C2:30
## Mode :character D3: 2
## D4: 1
## D5: 9
##
##
## SL_CM FL_CM TL_CM HH_CM
## Min. :16.90 Min. :14.00 Min. :14.90 Min. :3.400
## 1st Qu.:18.70 1st Qu.:16.50 1st Qu.:17.90 1st Qu.:3.700
## Median :19.40 Median :19.90 Median :21.60 Median :4.000
## Mean :20.98 Mean :20.76 Mean :22.81 Mean :4.245
## 3rd Qu.:24.00 3rd Qu.:25.20 3rd Qu.:28.00 3rd Qu.:4.800
## Max. :27.00 Max. :29.90 Max. :34.50 Max. :5.600
## NA's :48 NA's :50
## ED_CM W_G finclip otoliths gutted
## Min. :1.100 Min. : 52.7 Length:61 Length:61 N:61
## 1st Qu.:1.500 1st Qu.: 81.4 Class :character Class :character
## Median :1.600 Median :132.2 Mode :character Mode :character
## Mean :1.618 Mean :182.6
## 3rd Qu.:1.800 3rd Qu.:281.2
## Max. :2.100 Max. :481.6
## NA's :50
## gonads_present sex maturity G1L_CM G2L_CM
## Length:61 :12 :49 Min. :1.000 Min. :0.500
## Class :character F:33 B: 4 1st Qu.:2.500 1st Qu.:2.300
## Mode :character M:16 C: 8 Median :3.500 Median :3.200
## Mean :3.668 Mean :3.532
## 3rd Qu.:4.600 3rd Qu.:4.800
## Max. :6.500 Max. :6.500
## NA's :8 NA's :8
## GW_G stomach._sample
## Length:61 Length:61
## Class :character Class :character
## Mode :character Mode :character
##
##
##
## First look at the data with a scatterplot
fill = sex
ggplot(vsdf_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 = "Rhomboplites aurorubens, n = 61") +
scale_fill_viridis(option = "D", discrete = TRUE)
LM Linear model FL-W
vs2506_FLW_lm <- lm(FL_CM ~ W_G, data = vsdf_2506_FL_W)
summary(vs2506_FLW_lm)##
## Call:
## lm(formula = FL_CM ~ W_G, data = vsdf_2506_FL_W)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.95025 -0.61850 -0.07886 0.43363 2.43369
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.499954 0.258763 52.17 <2e-16 ***
## W_G 0.039764 0.001194 33.32 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.09 on 59 degrees of freedom
## Multiple R-squared: 0.9495, Adjusted R-squared: 0.9487
## F-statistic: 1110 on 1 and 59 DF, p-value: < 2.2e-16Scatterplot the residuals against the fitted values to check homogeneity of variance
plot(resid(vs2506_FLW_lm) ~ fitted(vs2506_FLW_lm))
The scatter plot looks like it indicates a nonlinear relationship
Plot a histogram to check normality of residuals
hist(resid(vs2506_FLW_lm))
Plot the model
vs2506_FLW_lm_plottemp <- ggplot(vsdf_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 = "Rhomboplites aurorubens, 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
vs2506_FLW_lm_plot <- vs2506_FLW_lm_plottemp +
geom_text(x = 14, y = 450,
label = expression(W~'(g)' == "0.040" * FL~'(cm)' + "13.500"~"," ~ r^{2} == 0.937 ~ "(Linear)"),
colour = safe_colorblind_palette[1],
hjust = 0)
vs2506_FLW_lm_plot
NLS Models
NLS 1
vs2506_FLW_nls <- nls_multstart(W_G~A*FL_CM^B,
data = vsdf_2506_FL_W,
iter = 500,
start_lower = c(A=0, B=2),
start_upper = c(A=0.5, B=3.5))
summary(vs2506_FLW_nls)##
## Formula: W_G ~ A * FL_CM^B
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## A 0.025612 0.005767 4.441 3.99e-05 ***
## B 2.879533 0.069435 41.471 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 17.49 on 59 degrees of freedom
##
## Number of iterations to convergence: 18
## Achieved convergence tolerance: 1.49e-08Plot a histogram of the residuals to check normality of residuals
#check residuals
hist(resid(vs2506_FLW_nls))
Calculate R-squared
rss_fitted <- sum(residuals(vs2506_FLW_nls)^2)
rss_constant <- sum((vsdf_2506_FL_W$W_G - mean(vsdf_2506_FL_W$W_G))^2)
approx_r_squared <- 1 - rss_fitted / rss_constant
approx_r_squared## [1] 0.9783608Calculate AIC
AIC(vs2506_FLW_nls)## [1] 526.2267Calculate the predicted W_kg values based on this NLS model
vs2506_FLW_nls_x <- c(14:30) #rounded actual min and max FL_CM values
vs2506_FLW_nls_predic <- data.frame(vs2506_FLW_nls_x)
vs2506_FLW_nls_f <- function(vs2506_FLW_nls_x) {
#formula using calculated parameters
return(0.025612*vs2506_FLW_nls_x^2.879533)
}
vs2506_FLW_nls_predic$vs2506_FLW_nls_y <- sapply(vs2506_FLW_nls_predic$vs2506_FLW_nls_x, vs2506_FLW_nls_f)
View(vs2506_FLW_nls_predic)Plot the model on the same graph as the other model
vs2506_FLW_nls_plottemp <- vs2506_FLW_lm_plot +
geom_smooth(aes(x = vs2506_FLW_nls_x, y = vs2506_FLW_nls_y),
data = vs2506_FLW_nls_predic,
method = "loess",
se = F,
color = safe_colorblind_palette[2],
linewidth = 0.75,
alpha = 0.75)
# add label with equation
vs2506_FLW_nls_plot <- vs2506_FLW_nls_plottemp +
geom_text(x = 14, y = 420,
label = expression(W~'(g)' == 2.561%*%10^{-02} * FL~'(cm)'^{"2.880"}~"," ~ R^{2} == 0.978 ~ '(NLS)'),
colour = safe_colorblind_palette[2],
hjust = 0)
vs2506_FLW_nls_plot
Brms1: with nls parameters and a wide normal distribution around them as priors
set priors and run model
#set priors
vs2506_FLW_prior1 <- prior(normal(0.025612, 0.005767) ,nlpar = "A")+prior(normal(2.879533, 0.069435),nlpar = "B")
vs2506_FLW_brm1 <- brm(bf(W_G~A*FL_CM^B, A + B ~ 1, nl = TRUE),data = vsdf_2506_FL_W,
prior = vs2506_FLW_prior1,
warmup = 1000,
iter = 5000,
chains = 4,
cores = 4,
control = list(adapt_delta = 0.8))
print(summary(vs2506_FLW_brm1), digits = 6)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: W_G ~ A * FL_CM^B
## A ~ 1
## B ~ 1
## Data: vsdf_2506_FL_W (Number of observations: 61)
## 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.026056 0.003402 0.019871 0.033137 1.002427 4435 4932
## B_Intercept 2.876753 0.040404 2.799834 2.957858 1.002290 4428 4965
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 17.801317 1.654917 14.917995 21.415191 1.000458 6533 5423
##
## 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(vs2506_FLW_brm1))
Calculate R-squared
rss_fitted <- sum(residuals(vs2506_FLW_brm1)^2)
rss_constant <- sum((vsdf_2506_FL_W$W_G - mean(vsdf_2506_FL_W$W_G))^2)
approx_r_squared <- 1 - rss_fitted / rss_constant
approx_r_squared ## [1] 0.726672smaller than the nls r-squared
Calculate the predicted W_G values based on this NLS model
vs2506_FLW_brm1_x <- c(14:30) #rounded actual min and max FL values
vs2506_FLW_brm1_predic <- data.frame(vs2506_FLW_brm1_x)
vs2506_FLW_brm1_f <- function(vs2506_FLW_brm1_x) {
#formula using calculated parameters
return(0.025989*vs2506_FLW_brm1_x^2.877509)
}
vs2506_FLW_brm1_predic$vs2506_FLW_brm1_y <- sapply(vs2506_FLW_brm1_predic$vs2506_FLW_brm1_x, vs2506_FLW_brm1_f)
View(vs2506_FLW_brm1_predic)Plot
vs2506_FLW_brm1_plottemp <- vs2506_FLW_nls_plot +
geom_smooth(aes(x = vs2506_FLW_brm1_x, y = vs2506_FLW_brm1_y),
data = vs2506_FLW_brm1_predic,
method = "loess",
se = F,
color = safe_colorblind_palette[4],
linewidth = 0.75,
alpha = 0.75)
# add label with equation
vs2506_FLW_brm1_plot <- vs2506_FLW_brm1_plottemp +
geom_text(x = 14, y = 390,
label = expression(W~'(g)' == 2.599%*%10^{-02} * FL~'(cm)'^{2.878}~"," ~ R^{2} == 0.726 ~ '(BRMS)'),
colour = safe_colorblind_palette[4],
hjust = 0)
vs2506_FLW_brm1_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
vs2506_FLW_nls_plot_title <- expression(italic("R. aurorubens")~"FL-W Relationship (n = 61)")
vs2506_FLW_nls_plottemp <- ggplot(vsdf_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 = vs2506_FLW_nls_plot_title) +
theme_classic() +
geom_smooth(aes(x = vs2506_FLW_nls_x, y = vs2506_FLW_nls_y),
data = vs2506_FLW_nls_predic,
method = "loess",
se = F,
color = turbo_colors_onboard[9],
linewidth = 0.75,
alpha = 0.75) +
ylim(0, 1600) +
xlim(10, 50)
# add label with equation
vs2506_FLW_nls_plot <- vs2506_FLW_nls_plottemp +
geom_label(x = 16, y = 1500,
label = expression(W~'(g)' == 2.561%*%10^{-02} * FL~'(cm)'^{"2.880"}~"," ~ R^{2} == 0.839 ~ '(NLS)'),
fill = "#ffa98e",
hjust = 0)
vs2506_FLW_nls_plot
Save plot
ggsave("vs2506_FLW_nls_plot.png", plot = vs2506_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
vsdf_2506_TL_W <- read.csv("vs2506df.csv")
vsdf_2506_TL_W$catch_month <- factor(vsdf_2506_TL_W$catch_month)
vsdf_2506_TL_W$year <- factor(vsdf_2506_TL_W$year)
vsdf_2506_TL_W$species <- factor(vsdf_2506_TL_W$species)
vsdf_2506_TL_W$common_name <- factor(vsdf_2506_TL_W$common_name)
vsdf_2506_TL_W$frozen <- factor(vsdf_2506_TL_W$frozen)
vsdf_2506_TL_W$catch_method <- factor(vsdf_2506_TL_W$catch_method)
vsdf_2506_TL_W$bait <- factor(vsdf_2506_TL_W$bait)
vsdf_2506_TL_W$area <- factor(vsdf_2506_TL_W$area)
vsdf_2506_TL_W$gutted <- factor(vsdf_2506_TL_W$gutted)
vsdf_2506_TL_W$sex <- factor(vsdf_2506_TL_W$sex)
vsdf_2506_TL_W$maturity <- factor(vsdf_2506_TL_W$maturity)
# Remove gutted samples
vsdf_2506_TL_W <- subset(filter(vsdf_2506_TL_W, gutted == "N"))
# remove NA rows
vsdf_2506_TL_W <- vsdf_2506_TL_W%>%
drop_na(TL_CM)
View(vsdf_2506_TL_W)
summary(vsdf_2506_TL_W)## X catch_date catch_month landing_time
## Min. : 1 Length:61 Apr: 1 Length:61
## 1st Qu.:16 Class :character Jul: 9 Class :character
## Median :31 Mode :character Jun:15 Mode :character
## Mean :31 May: 5
## 3rd Qu.:46 Sep:31
## Max. :61
##
## dissection_date dissection_time year month
## Length:61 Length:61 2023:49 Length:61
## Class :character Class :character 2024:12 Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
## observer species common_name
## Length:61 Rhomboplites aurorubens:61 Vermilion Snapper:61
## Class :character
## Mode :character
##
##
##
##
## fish_code frozen catch_method bait area
## Length:61 N:12 :49 :49 :19
## Class :character Y:49 Redfish trap:12 Japanese bait:12 C2:30
## Mode :character D3: 2
## D4: 1
## D5: 9
##
##
## SL_CM FL_CM TL_CM HH_CM
## Min. :16.90 Min. :14.00 Min. :14.90 Min. :3.400
## 1st Qu.:18.70 1st Qu.:16.50 1st Qu.:17.90 1st Qu.:3.700
## Median :19.40 Median :19.90 Median :21.60 Median :4.000
## Mean :20.98 Mean :20.76 Mean :22.81 Mean :4.245
## 3rd Qu.:24.00 3rd Qu.:25.20 3rd Qu.:28.00 3rd Qu.:4.800
## Max. :27.00 Max. :29.90 Max. :34.50 Max. :5.600
## NA's :48 NA's :50
## ED_CM W_G finclip otoliths gutted
## Min. :1.100 Min. : 52.7 Length:61 Length:61 N:61
## 1st Qu.:1.500 1st Qu.: 81.4 Class :character Class :character
## Median :1.600 Median :132.2 Mode :character Mode :character
## Mean :1.618 Mean :182.6
## 3rd Qu.:1.800 3rd Qu.:281.2
## Max. :2.100 Max. :481.6
## NA's :50
## gonads_present sex maturity G1L_CM G2L_CM
## Length:61 :12 :49 Min. :1.000 Min. :0.500
## Class :character F:33 B: 4 1st Qu.:2.500 1st Qu.:2.300
## Mode :character M:16 C: 8 Median :3.500 Median :3.200
## Mean :3.668 Mean :3.532
## 3rd Qu.:4.600 3rd Qu.:4.800
## Max. :6.500 Max. :6.500
## NA's :8 NA's :8
## GW_G stomach._sample
## Length:61 Length:61
## Class :character Class :character
## Mode :character Mode :character
##
##
##
## First look at the data with a scatterplot
fill = sex
ggplot(vsdf_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 = "Rhomboplites aurorubens, n = 61") +
scale_fill_viridis(option = "D", discrete = TRUE)
Linear model
vs2506_TLW_lm1 <- lm(TL_CM ~ W_G, data = vsdf_2506_TL_W)
summary(vs2506_TLW_lm1)##
## Call:
## lm(formula = TL_CM ~ W_G, data = vsdf_2506_TL_W)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.2078 -0.9735 -0.2731 0.5579 6.3355
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.435084 0.331881 43.49 <2e-16 ***
## W_G 0.045866 0.001531 29.96 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.398 on 59 degrees of freedom
## Multiple R-squared: 0.9383, Adjusted R-squared: 0.9373
## F-statistic: 897.8 on 1 and 59 DF, p-value: < 2.2e-16Scatterplot the residuals against the fitted values to check homogeneity of variance
plot(resid(vs2506_TLW_lm1) ~ fitted(vs2506_TLW_lm1))
looks nonlinear
Plot a histogram to check normality of residuals
qqnorm(resid(vs2506_TLW_lm1))
qqline(resid(vs2506_TLW_lm1))
Plot the model
vs2506_TLW_lm1_plottemp <- ggplot(vsdf_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 = "R. aurorubens, n = 61") +
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
vs2506_TLW_lm1_plot <- vs2506_TLW_lm1_plottemp +
geom_text(x = 14, y = 480,
label = expression(W~'(g)' == 0.046 * TL~'(cm)' + 14.435~"," ~ r^{2} == 0.9118 ~ "(Linear)"),
colour = safe_colorblind_palette[1],
hjust = 0)
vs2506_TLW_lm1_plot
NLS
vs2506_TLW_nls <- nls_multstart(W_G~A*TL_CM^B,
data = vsdf_2506_TL_W,
iter = 500,
start_lower = c(A=0, B=2),
start_upper = c(A=0.5, B=3.5))
summary(vs2506_TLW_nls)##
## Formula: W_G ~ A * TL_CM^B
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## A 0.033884 0.009484 3.573 0.000711 ***
## B 2.705590 0.083646 32.346 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 22.38 on 59 degrees of freedom
##
## Number of iterations to convergence: 9
## Achieved convergence tolerance: 1.49e-08Plot a histogram of the residuals to check normality of residuals
qqnorm(resid(vs2506_TLW_nls))
qqline(resid(vs2506_TLW_nls))
Calculate R-squared
rss_fitted <- sum(residuals(vs2506_TLW_nls)^2)
rss_constant <- sum((vsdf_2506_TL_W$W_G - mean(vsdf_2506_TL_W$W_G))^2)
approx_r_squared <- 1 - rss_fitted / rss_constant
approx_r_squared## [1] 0.9645734Calculate AIC
AIC(vs2506_TLW_nls)## [1] 556.2969Calculate the predicted W_kg values based on this NLS model
vs2506_TLW_nls_x <- c(14:35) #rounded actual min and max TL_CM values
vs2506_TLW_nls_predic <- data.frame(vs2506_TLW_nls_x)
vs2506_TLW_nls_f <- function(vs2506_TLW_nls_x) {
#formula using calculated parameters
return(0.033884*vs2506_TLW_nls_x^2.705590)
}
vs2506_TLW_nls_predic$vs2506_TLW_nls_y <- sapply(vs2506_TLW_nls_predic$vs2506_TLW_nls_x, vs2506_TLW_nls_f)
View(vs2506_TLW_nls_predic)Plot the model on the same graph as the other model
vs2506_TLW_nls_plottemp <- vs2506_TLW_lm1_plot +
geom_smooth(aes(x = vs2506_TLW_nls_x, y = vs2506_TLW_nls_y),
data = vs2506_TLW_nls_predic,
method = "loess",
se = F,
color = safe_colorblind_palette[2],
linewidth = 0.75)
# add label with equation
vs2506_TLW_nls_plot <- vs2506_TLW_nls_plottemp +
geom_text(x = 14, y = 450,
label = expression(W~'(g)' == 3.3884%*%10^{-03} * TL~'(cm)'^{2.706}~"," ~ R^{2} == 0.965 ~ '(NLS)'),
colour = safe_colorblind_palette[2],
hjust = 0)
vs2506_TLW_nls_plot
Brms: with nls parameters to start
set priors and run model
#set priors
vs2506_TLW_prior1 <- prior(normal(0.033884, 0.009484) ,nlpar = "A")+prior(normal(2.705590, 0.083646),nlpar = "B")
vs2506_TLW_brm1 <- brm(bf(W_G~A*TL_CM^B, A + B ~ 1, nl = TRUE),data = vsdf_2506_TL_W,
prior = vs2506_TLW_prior1,
warmup = 1000,
iter = 5000,
chains = 4,
cores = 4,
control = list(adapt_delta = 0.8))
print(summary(vs2506_TLW_brm1), digits = 6)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: W_G ~ A * TL_CM^B
## A ~ 1
## B ~ 1
## Data: vsdf_2506_TL_W (Number of observations: 61)
## 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.034711 0.005552 0.024606 0.046088 1.000490 4920 5504
## B_Intercept 2.702118 0.048194 2.612793 2.801258 1.000451 4912 5499
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 22.747334 2.139713 19.024633 27.389647 1.001442 6678 5983
##
## 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(vs2506_TLW_brm1))
Calculate R-squared
rss_fitted <- sum(residuals(vs2506_TLW_brm1)^2)
rss_constant <- sum((vsdf_2506_TL_W$W_G - mean(vsdf_2506_TL_W$W_G))^2)
approx_r_squared <- 1 - rss_fitted / rss_constant
approx_r_squared ## [1] 0.5529828Calculate the predicted W_G values based on this NLS model
vs2506_TLW_brm1_x <- c(14:35) #rounded actual min and max TL values
vs2506_TLW_brm1_predic <- data.frame(vs2506_TLW_brm1_x)
vs2506_TLW_brm1_f <- function(vs2506_TLW_brm1_x) {
#formula using calculated parameters
return(0.034724*vs2506_TLW_brm1_x^2.701907)
}
vs2506_TLW_brm1_predic$vs2506_TLW_brm1_y <- sapply(vs2506_TLW_brm1_predic$vs2506_TLW_brm1_x, vs2506_TLW_brm1_f)
View(vs2506_TLW_brm1_predic)Plot
vs2506_TLW_brm1_plottemp <- vs2506_TLW_nls_plot +
geom_smooth(aes(x = vs2506_TLW_brm1_x, y = vs2506_TLW_brm1_y),
data = vs2506_TLW_brm1_predic,
method = "loess",
se = F,
color = safe_colorblind_palette[4],
linewidth = 0.75)
# add label with equation
vs2506_TLW_brm1_plot <- vs2506_TLW_brm1_plottemp +
geom_text(x = 14, y = 420,
label = expression(W~'(g)' == 3.472%*%10^{-03} * TL~'(cm)'^{2.702}~"," ~ R^{2} == 0.555 ~ '(BRMS)'),
colour = safe_colorblind_palette[4],
hjust = 0)
vs2506_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
vsdf_2506_SL_W <- read.csv("vs2506df.csv")
vsdf_2506_SL_W$catch_month <- factor(vsdf_2506_SL_W$catch_month)
vsdf_2506_SL_W$year <- factor(vsdf_2506_SL_W$year)
vsdf_2506_SL_W$species <- factor(vsdf_2506_SL_W$species)
vsdf_2506_SL_W$common_name <- factor(vsdf_2506_SL_W$common_name)
vsdf_2506_SL_W$frozen <- factor(vsdf_2506_SL_W$frozen)
vsdf_2506_SL_W$catch_method <- factor(vsdf_2506_SL_W$catch_method)
vsdf_2506_SL_W$bait <- factor(vsdf_2506_SL_W$bait)
vsdf_2506_SL_W$area <- factor(vsdf_2506_SL_W$area)
vsdf_2506_SL_W$gutted <- factor(vsdf_2506_SL_W$gutted)
vsdf_2506_SL_W$sex <- factor(vsdf_2506_SL_W$sex)
vsdf_2506_SL_W$maturity <- factor(vsdf_2506_SL_W$maturity)
# Remove gutted samples
vsdf_2506_SL_W <- subset(filter(vsdf_2506_SL_W, gutted == "N"))
# remove NA rows
vsdf_2506_SL_W <- vsdf_2506_SL_W%>%
drop_na(SL_CM)
View(vsdf_2506_SL_W)
summary(vsdf_2506_SL_W)## X catch_date catch_month landing_time
## Min. : 1.00 Length:13 Apr:1 Length:13
## 1st Qu.: 4.00 Class :character Jul:0 Class :character
## Median : 7.00 Mode :character Jun:6 Mode :character
## Mean :10.69 May:5
## 3rd Qu.:10.00 Sep:1
## Max. :61.00
##
## dissection_date dissection_time year month
## Length:13 Length:13 2023: 1 Length:13
## Class :character Class :character 2024:12 Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
## observer species common_name
## Length:13 Rhomboplites aurorubens:13 Vermilion Snapper:13
## Class :character
## Mode :character
##
##
##
##
## fish_code frozen catch_method bait area
## Length:13 N:12 : 1 : 1 :1
## Class :character Y: 1 Redfish trap:12 Japanese bait:12 C2:0
## Mode :character D3:2
## D4:1
## D5:9
##
##
## SL_CM FL_CM TL_CM HH_CM ED_CM
## Min. :16.90 Min. :19.80 Min. :21.6 Min. :3.400 Min. :1.100
## 1st Qu.:18.70 1st Qu.:21.30 1st Qu.:23.1 1st Qu.:3.700 1st Qu.:1.500
## Median :19.40 Median :21.70 Median :27.1 Median :4.000 Median :1.600
## Mean :20.98 Mean :23.99 Mean :26.7 Mean :4.245 Mean :1.618
## 3rd Qu.:24.00 3rd Qu.:27.00 3rd Qu.:29.6 3rd Qu.:4.800 3rd Qu.:1.800
## Max. :27.00 Max. :29.90 Max. :34.5 Max. :5.600 Max. :2.100
## NA's :2 NA's :2
## W_G finclip otoliths gutted
## Min. :131.0 Length:13 Length:13 N:13
## 1st Qu.:142.5 Class :character Class :character
## Median :170.6 Mode :character Mode :character
## Mean :241.7
## 3rd Qu.:328.6
## Max. :481.6
##
## gonads_present sex maturity G1L_CM G2L_CM
## Length:13 :0 :1 Min. :2.500 Min. :2.100
## Class :character F:8 B:4 1st Qu.:3.700 1st Qu.:3.200
## Mode :character M:5 C:8 Median :4.200 Median :4.500
## Mean :4.431 Mean :4.215
## 3rd Qu.:5.200 3rd Qu.:5.100
## Max. :6.500 Max. :5.900
##
## GW_G stomach._sample
## Length:13 Length:13
## Class :character Class :character
## Mode :character Mode :character
##
##
##
## First look at the data with a scatterplot
fill = sex
ggplot(vsdf_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 = "Rhomboplites aurorubens, n = 24") +
scale_fill_viridis(option = "D", discrete = TRUE)
LM Linear model SL-W
vs2506_FLW_lm <- lm(SL_CM ~ W_G, data = vsdf_2506_SL_W)
summary(vs2506_FLW_lm)##
## Call:
## lm(formula = SL_CM ~ W_G, data = vsdf_2506_SL_W)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.03817 -0.59112 0.09394 0.28339 1.17539
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.261730 0.401990 35.48 1.07e-12 ***
## W_G 0.027810 0.001488 18.69 1.10e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6475 on 11 degrees of freedom
## Multiple R-squared: 0.9695, Adjusted R-squared: 0.9667
## F-statistic: 349.4 on 1 and 11 DF, p-value: 1.102e-09Scatterplot the residuals against the fitted values to check homogeneity of variance
plot(resid(vs2506_FLW_lm) ~ fitted(vs2506_FLW_lm))
looks like a nonlinear relationship
Plot a histogram to check normality of residuals
qqnorm(resid(vs2506_FLW_lm))
qqline(resid(vs2506_FLW_lm))
not normal
Plot the model
vs2506_SLW_lm_plottemp <- ggplot(vsdf_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 = "Rhomboplites aurorubens, 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
vs2506_SLW_lm_plot <- vs2506_SLW_lm_plottemp +
geom_text(x = 16, y = 450,
label = expression(W~'(g)' == 0.028 * SL~'(cm)' + 14.262~"," ~ r^{2} == 0.9667 ~ "(Linear)"),
colour = safe_colorblind_palette[1],
hjust = 0)
vs2506_SLW_lm_plot
NLS Models
NLS 1
vs2506_SLW_nls <- nls_multstart(W_G~A*SL_CM^B,
data = vsdf_2506_SL_W,
iter = 500,
start_lower = c(A=0, B=2),
start_upper = c(A=0.5, B=3.5))
print(summary(vs2506_SLW_nls), digits = 3)##
## Formula: W_G ~ A * SL_CM^B
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## A 0.02408 0.00727 3.31 0.0069 **
## B 3.00100 0.09501 31.58 3.8e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 12.6 on 11 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(vs2506_SLW_nls))
Calculate R-squared
rss_fitted <- sum(residuals(vs2506_SLW_nls)^2)
rss_constant <- sum((vsdf_2506_SL_W$W_G - mean(vsdf_2506_SL_W$W_G))^2)
approx_r_squared <- 1 - rss_fitted / rss_constant
approx_r_squared## [1] 0.9907094Calculate AIC
AIC(vs2506_SLW_nls)## [1] 106.6995Calculate the predicted W_g values based on this NLS model
vs2506_SLW_nls_x <- c(16:27) #rounded actual min and max FL_CM values
vs2506_SLW_nls_predic <- data.frame(vs2506_SLW_nls_x)
vs2506_SLW_nls_f <- function(vs2506_SLW_nls_x) {
#formula using calculated parameters
return(0.02408*vs2506_SLW_nls_x^3.00100)
}
vs2506_SLW_nls_predic$vs2506_SLW_nls_y <- sapply(vs2506_SLW_nls_predic$vs2506_SLW_nls_x, vs2506_SLW_nls_f)
View(vs2506_SLW_nls_predic)Plot the model on the same graph as the other model
vs2506_SLW_nls_plottemp <- vs2506_SLW_lm_plot +
geom_smooth(aes(x = vs2506_SLW_nls_x, y = vs2506_SLW_nls_y),
data = vs2506_SLW_nls_predic,
method = "loess",
se = F,
color = safe_colorblind_palette[2],
linewidth = 0.75,
alpha = 0.75)
# add label with equation
vs2506_SLW_nls_plot <- vs2506_SLW_nls_plottemp +
geom_text(x = 16, y = 420,
label = expression(W~'(g)' == 2.408%*%10^{-02} * SL~'(cm)'^{3.001}~"," ~ R^{2} == 0.991 ~ '(NLS)'),
colour = safe_colorblind_palette[2],
hjust = 0)
vs2506_SLW_nls_plot
Brms: with nls parameters to start
set priors and run model
#set priors
vs2506_SLW_prior1 <- prior(normal(0.02408, 0.00727) ,nlpar = "A")+prior(normal(3.00100, 0.09501),nlpar = "B")
vs2506_SLW_brm1 <- brm(bf(W_G~A*SL_CM^B, A + B ~ 1, nl = TRUE),data = vsdf_2506_SL_W,
prior = vs2506_SLW_prior1,
warmup = 1000,
iter = 5000,
chains = 4,
cores = 4,
control = list(adapt_delta = 0.8))
print(summary(vs2506_SLW_brm1), digits = 6)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: W_G ~ A * SL_CM^B
## A ~ 1
## B ~ 1
## Data: vsdf_2506_SL_W (Number of observations: 13)
## 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.024798 0.004286 0.017161 0.033878 1.000839 3841 4246
## B_Intercept 2.996383 0.054696 2.893468 3.108037 1.000847 3843 4220
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 13.714448 3.196995 9.075373 21.382235 1.000366 5104 5375
##
## 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(vs2506_SLW_brm1))
Calculate R-squared
rss_fitted <- sum(residuals(vs2506_SLW_brm1)^2)
rss_constant <- sum((vsdf_2506_SL_W$W_G - mean(vsdf_2506_SL_W$W_G))^2)
approx_r_squared <- 1 - rss_fitted / rss_constant
approx_r_squared ## [1] 0.8388968Calculate the predicted W_G values based on this NLS model
vs2506_SLW_brm1_x <- c(16:27) #rounded actual min and max SL values
vs2506_SLW_brm1_predic <- data.frame(vs2506_SLW_brm1_x)
vs2506_SLW_brm1_f <- function(vs2506_SLW_brm1_x) {
#formula using calculated parameters
return(0.024737*vs2506_SLW_brm1_x^2.996963)
}
vs2506_SLW_brm1_predic$vs2506_SLW_brm1_y <- sapply(vs2506_SLW_brm1_predic$vs2506_SLW_brm1_x, vs2506_SLW_brm1_f)
View(vs2506_SLW_brm1_predic)Plot
vs2506_SLW_brm1_plottemp <- vs2506_SLW_nls_plot +
geom_smooth(aes(x = vs2506_SLW_brm1_x, y = vs2506_SLW_brm1_y),
data = vs2506_SLW_brm1_predic,
method = "loess",
se = F,
color = safe_colorblind_palette[3],
linewidth = 0.75)
# add label with equation
vs2506_SLW_brm1_plot <- vs2506_SLW_brm1_plottemp +
geom_text(x = 16, y = 390,
label = expression(W~'(g)' == 2.474%*%10^{-02} * SL~'(cm)'^{2.997}~"," ~ R^{2} == 0.839 ~ '(BRMS)'),
colour = safe_colorblind_palette[3],
hjust = 0)
vs2506_SLW_brm1_plot
Final decision: NLS
Maturity
Load data file
vsdf_2506_maturity <- read.csv("vs2506df.csv")
vsdf_2506_maturity$catch_month <- factor(vsdf_2506_maturity$catch_month)
vsdf_2506_maturity$year <- factor(vsdf_2506_maturity$year)
vsdf_2506_maturity$species <- factor(vsdf_2506_maturity$species)
vsdf_2506_maturity$common_name <- factor(vsdf_2506_maturity$common_name)
vsdf_2506_maturity$frozen <- factor(vsdf_2506_maturity$frozen)
vsdf_2506_maturity$catch_method <- factor(vsdf_2506_maturity$catch_method)
vsdf_2506_maturity$bait <- factor(vsdf_2506_maturity$bait)
vsdf_2506_maturity$area <- factor(vsdf_2506_maturity$area)
vsdf_2506_maturity$gutted <- factor(vsdf_2506_maturity$gutted)
vsdf_2506_maturity$sex <- factor(vsdf_2506_maturity$sex)
vsdf_2506_maturity$maturity <- factor(vsdf_2506_maturity$maturity)
#drop the columns with neither sex nor maturity measurements
vsdf_2506_maturity <- vsdf_2506_maturity %>%
filter(maturity != "" | sex != "")calculate sex ratio
vs_totalfemale <- vsdf_2506_maturity %>%
filter(sex == "F") %>%
count()
vs_totalmale <- vsdf_2506_maturity %>%
filter(sex == "M") %>%
count()
vs_sexratio <- vs_totalfemale/vs_totalmale
vs_sexratio #female number## n
## 1 2.0625Descriptive statistics
vs_sex_summary <- vsdf_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)
)
vs_sex_summary## # A tibble: 2 × 8
## sex count mean sd median min max IQR
## <fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 F 33 20.8 4.87 19.9 14.5 29.9 8.3
## 2 M 16 20.7 4.32 20.0 15.2 28 7.85vsdf_2506_maturity <- vsdf_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))
vs_maturity_summary <- vsdf_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)
)
vs_maturity_summary## # A tibble: 2 × 8
## maturity_state count mean sd median min max IQR
## <fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 mature 12 23.7 3.70 21.6 19.8 29.9 4.65
## 2 <NA> 37 19.8 4.55 17.5 14.5 29.4 7.1