L. synagris body size models
Naomi Bakken
2024-08-15
This R markdown file includes all the R code used to analyze Lane Snapper (Lutjanus sayngris) 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)
library(MASS) #box cox lambda estimationsDefine 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
lsdf_2506_LLrelationships <- read.csv("ls2506df.csv")
lsdf_2506_LLrelationships$catch_month <- factor(lsdf_2506_LLrelationships$catch_month)
lsdf_2506_LLrelationships$year <- factor(lsdf_2506_LLrelationships$year)
lsdf_2506_LLrelationships$species <- factor(lsdf_2506_LLrelationships$species)
lsdf_2506_LLrelationships$common_name <- factor(lsdf_2506_LLrelationships$common_name)
lsdf_2506_LLrelationships$frozen <- factor(lsdf_2506_LLrelationships$frozen)
lsdf_2506_LLrelationships$catch_method <- factor(lsdf_2506_LLrelationships$catch_method)
lsdf_2506_LLrelationships$bait <- factor(lsdf_2506_LLrelationships$bait)
lsdf_2506_LLrelationships$area <- factor(lsdf_2506_LLrelationships$area)
lsdf_2506_LLrelationships$gutted <- factor(lsdf_2506_LLrelationships$gutted)
lsdf_2506_LLrelationships$sex <- factor(lsdf_2506_LLrelationships$sex)
lsdf_2506_LLrelationships$maturity <- factor(lsdf_2506_LLrelationships$maturity)
summary(lsdf_2506_LLrelationships)## X catch_date catch_month landing_time
## Min. : 1.00 Length:30 Apr: 1 Length:30
## 1st Qu.: 8.25 Class :character Aug: 3 Class :character
## Median :15.50 Mode :character Jul: 7 Mode :character
## Mean :15.50 Jun: 4
## 3rd Qu.:22.75 May:15
## Max. :30.00
##
## dissection_date dissection_time year month
## Length:30 Length:30 2023:11 Length:30
## Class :character Class :character 2024:19 Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
## observer species common_name fish_code
## Length:30 Lutjanus synagris:30 Lane Snapper:30 Length:30
## Class :character Class :character
## Mode :character Mode :character
##
##
##
##
## frozen catch_method bait area SL_CM
## N:19 :11 :11 :11 Min. :13.70
## Y:11 Lobster trap:19 Cowskin:19 B3 : 1 1st Qu.:19.00
## B3/4 : 2 Median :23.40
## B4 : 9 Mean :22.55
## B5 : 2 3rd Qu.:26.60
## D3 : 2 Max. :28.50
## LEAYLE: 3 NA's :11
## FL_CM TL_CM HH_CM ED_CM
## Min. :16.10 Min. :16.90 Min. :2.900 Min. :0.800
## 1st Qu.:23.82 1st Qu.:25.60 1st Qu.:4.200 1st Qu.:1.100
## Median :26.65 Median :28.75 Median :5.500 Median :1.250
## Mean :26.28 Mean :28.05 Mean :5.292 Mean :1.283
## 3rd Qu.:30.35 3rd Qu.:32.40 3rd Qu.:6.175 3rd Qu.:1.550
## Max. :33.50 Max. :35.10 Max. :7.000 Max. :1.900
## NA's :12 NA's :12
## W_G finclip otoliths gutted
## Min. : 65.7 Length:30 Length:30 N:30
## 1st Qu.:212.2 Class :character Class :character
## Median :315.2 Mode :character Mode :character
## Mean :319.8
## 3rd Qu.:454.9
## Max. :563.2
##
## gonads_present sex maturity G1L_CM G2L_CM
## Length:30 F:20 :11 Min. :2.000 Min. :1.200
## Class :character M:10 A: 1 1st Qu.:3.825 1st Qu.:4.000
## Mode :character B: 3 Median :4.750 Median :4.900
## C:15 Mean :4.630 Mean :4.910
## 3rd Qu.:5.375 3rd Qu.:5.975
## Max. :7.700 Max. :8.200
##
## GW_G stomach._sample
## Length:30 Length:30
## Class :character Class :character
## Mode :character Mode :character
##
##
##
## View(lsdf_2506_LLrelationships)
sumtable(lsdf_2506_LLrelationships, digits = 4)| Variable | N | Mean | Std. Dev. | Min | Pctl. 25 | Pctl. 75 | Max |
|---|---|---|---|---|---|---|---|
| X | 30 | 15.5 | 8.803 | 1 | 8.25 | 22.75 | 30 |
| catch_month | 30 | ||||||
| … Apr | 1 | 3.33% | |||||
| … Aug | 3 | 10% | |||||
| … Jul | 7 | 23.33% | |||||
| … Jun | 4 | 13.33% | |||||
| … May | 15 | 50% | |||||
| year | 30 | ||||||
| … 2023 | 11 | 36.67% | |||||
| … 2024 | 19 | 63.33% | |||||
| month | 30 | ||||||
| … Apr | 1 | 3.33% | |||||
| … Aug | 3 | 10% | |||||
| … Jul | 7 | 23.33% | |||||
| … Jun | 4 | 13.33% | |||||
| … May | 15 | 50% | |||||
| observer | 30 | ||||||
| … HK | 7 | 23.33% | |||||
| … NB | 1 | 3.33% | |||||
| … QH | 4 | 13.33% | |||||
| … SH | 4 | 13.33% | |||||
| … SM | 2 | 6.67% | |||||
| … TO | 12 | 40% | |||||
| species | 30 | ||||||
| … Lutjanus synagris | 30 | 100% | |||||
| common_name | 30 | ||||||
| … Lane Snapper | 30 | 100% | |||||
| frozen | 30 | ||||||
| … N | 19 | 63.33% | |||||
| … Y | 11 | 36.67% | |||||
| catch_method | 30 | ||||||
| … | 11 | 36.67% | |||||
| … Lobster trap | 19 | 63.33% | |||||
| bait | 30 | ||||||
| … | 11 | 36.67% | |||||
| … Cowskin | 19 | 63.33% | |||||
| area | 30 | ||||||
| … | 11 | 36.67% | |||||
| … B3 | 1 | 3.33% | |||||
| … B3/4 | 2 | 6.67% | |||||
| … B4 | 9 | 30% | |||||
| … B5 | 2 | 6.67% | |||||
| … D3 | 2 | 6.67% | |||||
| … LEAYLE | 3 | 10% | |||||
| SL_CM | 19 | 22.55 | 4.884 | 13.7 | 19 | 26.6 | 28.5 |
| FL_CM | 30 | 26.28 | 4.774 | 16.1 | 23.82 | 30.35 | 33.5 |
| TL_CM | 30 | 28.05 | 5.176 | 16.9 | 25.6 | 32.4 | 35.1 |
| HH_CM | 18 | 5.292 | 1.191 | 2.9 | 4.2 | 6.175 | 7 |
| ED_CM | 18 | 1.283 | 0.3167 | 0.8 | 1.1 | 1.55 | 1.9 |
| W_G | 30 | 319.8 | 151.8 | 65.7 | 212.2 | 454.9 | 563.2 |
| finclip | 30 | ||||||
| … | 11 | 36.67% | |||||
| … N | 2 | 6.67% | |||||
| … Y | 17 | 56.67% | |||||
| otoliths | 30 | ||||||
| … Y | 30 | 100% | |||||
| gutted | 30 | ||||||
| … N | 30 | 100% | |||||
| gonads_present | 30 | ||||||
| … | 1 | 3.33% | |||||
| … y | 1 | 3.33% | |||||
| … Y | 28 | 93.33% | |||||
| sex | 30 | ||||||
| … F | 20 | 66.67% | |||||
| … M | 10 | 33.33% | |||||
| maturity | 30 | ||||||
| … | 11 | 36.67% | |||||
| … A | 1 | 3.33% | |||||
| … B | 3 | 10% | |||||
| … C | 15 | 50% | |||||
| G1L_CM | 30 | 4.63 | 1.358 | 2 | 3.825 | 5.375 | 7.7 |
| G2L_CM | 30 | 4.91 | 1.585 | 1.2 | 4 | 5.975 | 8.2 |
| stomach._sample | 30 | ||||||
| … N | 29 | 96.67% | |||||
| … Y | 1 | 3.33% |
Fork length (FL_CM) - Total length (TL_CM) relationship
ls2506_FLTL_plotbase <- ggplot(lsdf_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 synagris FL-TL relationship (n = 82)",
fill = "Sex") +
scale_fill_viridis(option = "C",
discrete = TRUE,
begin = 0.2,
end = 0.8,
labels = c("N/A", "Female", "Male"))
ls2506_FLTL_plotbase
Linear model
ls2506_FLTL_lm <- lm(FL_CM ~ TL_CM, data = lsdf_2506_LLrelationships)
summary(ls2506_FLTL_lm)##
## Call:
## lm(formula = FL_CM ~ TL_CM, data = lsdf_2506_LLrelationships)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.47691 -0.26079 -0.02869 0.18302 0.82095
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.46908 0.33414 1.404 0.171
## TL_CM 0.92028 0.01172 78.508 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3267 on 28 degrees of freedom
## Multiple R-squared: 0.9955, Adjusted R-squared: 0.9953
## F-statistic: 6164 on 1 and 28 DF, p-value: < 2.2e-16length(ls2506_FLTL_lm$fitted.values)## [1] 30Scatterplot the residuals against the fitted values to check homogeneity of variance
plot(resid(ls2506_FLTL_lm) ~ fitted(ls2506_FLTL_lm))
Fine
Plot a histogram to check normality of residuals
qqnorm(ls2506_FLTL_lm$residuals)
qqline(ls2506_FLTL_lm$residuals)
slight positive skew - try a data transformation
Transform data
sqrtFL_CM <- sqrt(lsdf_2506_LLrelationships$FL_CM)
lsdf_2506_LLrelationships$lnFL_CM <- sqrtFL_CM
sqrtTL_CM <- sqrt(lsdf_2506_LLrelationships$TL_CM)
lsdf_2506_LLrelationships$lnTL_CM <- sqrtTL_CM
lnFL_CM <- log(lsdf_2506_LLrelationships$FL_CM)
lsdf_2506_LLrelationships$lnFL_CM <- lnFL_CM
lnTL_CM <- log(lsdf_2506_LLrelationships$TL_CM)
lsdf_2506_LLrelationships$lnTL_CM <- lnTL_CMLinear model sqrt transformation
ls2506_sqrt_FLTL_lm <- lm(sqrtFL_CM ~ sqrtTL_CM, data = lsdf_2506_LLrelationships)
summary(ls2506_sqrt_FLTL_lm)##
## Call:
## lm(formula = sqrtFL_CM ~ sqrtTL_CM, data = lsdf_2506_LLrelationships)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.04781 -0.02562 -0.00176 0.01933 0.07413
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.11436 0.06106 1.873 0.0716 .
## sqrtTL_CM 0.94648 0.01153 82.090 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03155 on 28 degrees of freedom
## Multiple R-squared: 0.9959, Adjusted R-squared: 0.9957
## F-statistic: 6739 on 1 and 28 DF, p-value: < 2.2e-16length(ls2506_sqrt_FLTL_lm$fitted.values)## [1] 30Scatterplot the residuals against the fitted values to check homogeneity of variance
plot(resid(ls2506_sqrt_FLTL_lm) ~ fitted(ls2506_sqrt_FLTL_lm))
Fine
Plot a histogram to check normality of residuals
qqnorm(resid(ls2506_sqrt_FLTL_lm))
qqline(resid(ls2506_sqrt_FLTL_lm))
Fork length (FL_CM) - Standard length (SL_CM) relationship
ls2506_FLSL_plot <- ggplot(lsdf_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 synagris FL-SL relationship (n = 19)",
fill = "Sex") +
scale_fill_viridis(option = "C",
discrete = TRUE,
begin = 0.2,
end = 0.8,
labels = c("N/A", "Female", "Male"))
ls2506_FLSL_plot
Linear model
ls2506_FLSL_lm <- lm(FL_CM ~ SL_CM, data = lsdf_2506_LLrelationships)
summary(ls2506_FLSL_lm)##
## Call:
## lm(formula = FL_CM ~ SL_CM, data = lsdf_2506_LLrelationships)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.5208 -0.1795 -0.0383 0.2652 0.5572
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.97341 0.35969 2.706 0.015 *
## SL_CM 1.10295 0.01561 70.660 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3234 on 17 degrees of freedom
## (11 observations deleted due to missingness)
## Multiple R-squared: 0.9966, Adjusted R-squared: 0.9964
## F-statistic: 4993 on 1 and 17 DF, p-value: < 2.2e-16length(ls2506_FLSL_lm$fitted.values)## [1] 19Scatterplot the residuals against the fitted values to check homogeneity of variance
plot(resid(ls2506_FLSL_lm) ~ fitted(ls2506_FLSL_lm))
Plot a histogram to check normality of residuals
qqnorm(resid(ls2506_FLSL_lm))
qqline(resid(ls2506_FLSL_lm))
Total length (FL_CM) - Standard length (SL_CM) relationship
ls2506_TLSL_plot <- ggplot(lsdf_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 synagris TL-SL relationship (n = 19)",
fill = "Sex") +
scale_fill_viridis(option = "C",
discrete = TRUE,
begin = 0.2,
end = 0.8,
labels = c("N/A", "Female", "Male"))
ls2506_FLSL_plot
Linear model
ls2506_TLSL_lm <- lm(TL_CM ~ SL_CM, data = lsdf_2506_LLrelationships)
summary(ls2506_TLSL_lm)##
## Call:
## lm(formula = TL_CM ~ SL_CM, data = lsdf_2506_LLrelationships)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.86831 -0.18799 0.00189 0.10724 0.91574
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.42140 0.43059 0.979 0.341
## SL_CM 1.20213 0.01869 64.332 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3872 on 17 degrees of freedom
## (11 observations deleted due to missingness)
## Multiple R-squared: 0.9959, Adjusted R-squared: 0.9957
## F-statistic: 4139 on 1 and 17 DF, p-value: < 2.2e-16length(ls2506_TLSL_lm$fitted.values)## [1] 19Scatterplot the residuals against the fitted values to check homogeneity of variance
plot(resid(ls2506_TLSL_lm) ~ fitted(ls2506_TLSL_lm))
Plot a histogram to check normality of residuals
qqnorm(resid(ls2506_TLSL_lm))
qqline(resid(ls2506_TLSL_lm))
heavy tails but fine
Fork length (FL_CM) - Weight (W_G) relationship
Load data file
lsdf_2506_FL_W <- read.csv("ls2506df.csv")
lsdf_2506_FL_W$catch_month <- factor(lsdf_2506_FL_W$catch_month)
lsdf_2506_FL_W$year <- factor(lsdf_2506_FL_W$year)
lsdf_2506_FL_W$species <- factor(lsdf_2506_FL_W$species)
lsdf_2506_FL_W$common_name <- factor(lsdf_2506_FL_W$common_name)
lsdf_2506_FL_W$frozen <- factor(lsdf_2506_FL_W$frozen)
lsdf_2506_FL_W$catch_method <- factor(lsdf_2506_FL_W$catch_method)
lsdf_2506_FL_W$bait <- factor(lsdf_2506_FL_W$bait)
lsdf_2506_FL_W$area <- factor(lsdf_2506_FL_W$area)
lsdf_2506_FL_W$gutted <- factor(lsdf_2506_FL_W$gutted)
lsdf_2506_FL_W$sex <- factor(lsdf_2506_FL_W$sex)
lsdf_2506_FL_W$maturity <- factor(lsdf_2506_FL_W$maturity)
# Remove gutted samples
lsdf_2506_FL_W <- subset(filter(lsdf_2506_FL_W, gutted == "N"))
# remove NA rows
lsdf_2506_FL_W <- lsdf_2506_FL_W%>%
drop_na(FL_CM)
View(lsdf_2506_FL_W)
summary(lsdf_2506_FL_W)## X catch_date catch_month landing_time
## Min. : 1.00 Length:30 Apr: 1 Length:30
## 1st Qu.: 8.25 Class :character Aug: 3 Class :character
## Median :15.50 Mode :character Jul: 7 Mode :character
## Mean :15.50 Jun: 4
## 3rd Qu.:22.75 May:15
## Max. :30.00
##
## dissection_date dissection_time year month
## Length:30 Length:30 2023:11 Length:30
## Class :character Class :character 2024:19 Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
## observer species common_name fish_code
## Length:30 Lutjanus synagris:30 Lane Snapper:30 Length:30
## Class :character Class :character
## Mode :character Mode :character
##
##
##
##
## frozen catch_method bait area SL_CM
## N:19 :11 :11 :11 Min. :13.70
## Y:11 Lobster trap:19 Cowskin:19 B3 : 1 1st Qu.:19.00
## B3/4 : 2 Median :23.40
## B4 : 9 Mean :22.55
## B5 : 2 3rd Qu.:26.60
## D3 : 2 Max. :28.50
## LEAYLE: 3 NA's :11
## FL_CM TL_CM HH_CM ED_CM
## Min. :16.10 Min. :16.90 Min. :2.900 Min. :0.800
## 1st Qu.:23.82 1st Qu.:25.60 1st Qu.:4.200 1st Qu.:1.100
## Median :26.65 Median :28.75 Median :5.500 Median :1.250
## Mean :26.28 Mean :28.05 Mean :5.292 Mean :1.283
## 3rd Qu.:30.35 3rd Qu.:32.40 3rd Qu.:6.175 3rd Qu.:1.550
## Max. :33.50 Max. :35.10 Max. :7.000 Max. :1.900
## NA's :12 NA's :12
## W_G finclip otoliths gutted
## Min. : 65.7 Length:30 Length:30 N:30
## 1st Qu.:212.2 Class :character Class :character
## Median :315.2 Mode :character Mode :character
## Mean :319.8
## 3rd Qu.:454.9
## Max. :563.2
##
## gonads_present sex maturity G1L_CM G2L_CM
## Length:30 F:20 :11 Min. :2.000 Min. :1.200
## Class :character M:10 A: 1 1st Qu.:3.825 1st Qu.:4.000
## Mode :character B: 3 Median :4.750 Median :4.900
## C:15 Mean :4.630 Mean :4.910
## 3rd Qu.:5.375 3rd Qu.:5.975
## Max. :7.700 Max. :8.200
##
## GW_G stomach._sample
## Length:30 Length:30
## Class :character Class :character
## Mode :character Mode :character
##
##
##
## First look at the data with a scatterplot
fill = sex
ggplot(lsdf_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 synagris, n = 30") +
scale_fill_viridis(option = "D", discrete = TRUE)
LM Linear model FL-W
ls2506_FLW_lm <- lm(FL_CM ~ W_G, data = lsdf_2506_FL_W)
summary(ls2506_FLW_lm)##
## Call:
## lm(formula = FL_CM ~ W_G, data = lsdf_2506_FL_W)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.3641 -0.8409 0.1387 0.7435 1.6429
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 16.443589 0.440903 37.30 <2e-16 ***
## W_G 0.030754 0.001249 24.62 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.021 on 28 degrees of freedom
## Multiple R-squared: 0.9558, Adjusted R-squared: 0.9543
## F-statistic: 606.1 on 1 and 28 DF, p-value: < 2.2e-16Scatterplot the residuals against the fitted values to check homogeneity of variance
plot(resid(ls2506_FLW_lm) ~ fitted(ls2506_FLW_lm))
The scatter plot clearly looks like it indicates a nonlinear
relationship
Plot a histogram to check normality of residuals
qqnorm(resid(ls2506_FLW_lm))
qqline(resid(ls2506_FLW_lm))
negative skew
Plot the model
ls2506_FLW_lm_plottemp <- ggplot(lsdf_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 synagris, n = 30") +
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
ls2506_FLW_lm_plot <- ls2506_FLW_lm_plottemp +
geom_text(x = 16, y = 550,
label = expression(W~'(g)' == 0.031 * TL~'(cm)' + 16.444~"," ~ R^{2} == 0.954 ~ "(Linear)"),
colour = safe_colorblind_palette[1],
hjust = 0)
ls2506_FLW_lm_plot
NLS Models
NLS 1 *
ls2506_FLW_nls <- nls_multstart(W_G~A*FL_CM^B,
data = lsdf_2506_FL_W,
iter = 500,
start_lower = c(A=0, B=2),
start_upper = c(A=0.5, B=3.5))
summary(ls2506_FLW_nls)##
## Formula: W_G ~ A * FL_CM^B
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## A 0.02487 0.00916 2.715 0.0112 *
## B 2.87101 0.10877 26.395 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 24.25 on 28 degrees of freedom
##
## Number of iterations to convergence: 24
## Achieved convergence tolerance: 1.49e-08Plot a histogram of the residuals to check normality of residuals
#check residuals
qqnorm(resid(ls2506_FLW_nls))
qqline(resid(ls2506_FLW_nls))
okay
Calculate R-squared
rss_fitted <- sum(residuals(ls2506_FLW_nls)^2)
rss_constant <- sum((lsdf_2506_FL_W$W_G - mean(lsdf_2506_FL_W$W_G))^2)
approx_r_squared <- 1 - rss_fitted / rss_constant
approx_r_squared## [1] 0.9753532Calculate AIC
AIC(ls2506_FLW_nls)## [1] 280.3699Calculate the predicted W_kg values based on this NLS model
ls2506_FLW_nls_x <- c(16:34) #rounded actual min and max FL_MM values
ls2506_FLW_nls_predic <- data.frame(ls2506_FLW_nls_x)
ls2506_FLW_nls_f <- function(ls2506_FLW_nls_x) {
#formula using calculated parameters
return(0.02487*ls2506_FLW_nls_x^2.87101)
}
ls2506_FLW_nls_predic$ls2506_FLW_nls_y <- sapply(ls2506_FLW_nls_predic$ls2506_FLW_nls_x, ls2506_FLW_nls_f)
View(ls2506_FLW_nls_predic)Plot the model on the same graph as the other model
ls2506_FLW_nls_plottemp <- ls2506_FLW_lm_plot +
geom_smooth(aes(x = ls2506_FLW_nls_x, y = ls2506_FLW_nls_y),
data = ls2506_FLW_nls_predic,
method = "loess",
se = F,
color = safe_colorblind_palette[2],
linewidth = 0.75,
alpha = 0.75)
# add label with equation
ls2506_FLW_nls_plot <- ls2506_FLW_nls_plottemp +
geom_text(x = 16, y = 510,
label = expression(W~'(g)' == 2.487%*%10^{-02} * FL~'(cm)'^{2.871}~"," ~ R^{2} == 0.975 ~ '(NLS)'),
colour = safe_colorblind_palette[2],
hjust = 0)
ls2506_FLW_nls_plot
Brms1: with nls parameters and a wide normal distribution around them as priors
set priors and run model
#set priors
ls2506_FLW_prior1 <- prior(normal(0.02487, 0.01) ,nlpar = "A")+prior(normal(2.87101, 1),nlpar = "B")
ls2506_FLW_brm1 <- brm(bf(W_G~A*FL_CM^B, A + B ~ 1, nl = TRUE),data = lsdf_2506_FL_W,
prior = ls2506_FLW_prior1,
warmup = 1000,
iter = 5000,
chains = 4,
cores = 4,
control = list(adapt_delta = 0.8))
print(summary(ls2506_FLW_brm1), digits = 6)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: W_G ~ A * FL_CM^B
## A ~ 1
## B ~ 1
## Data: lsdf_2506_FL_W (Number of observations: 30)
## 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.026130 0.006569 0.014623 0.040289 1.000413 3691 4383
## B_Intercept 2.865920 0.076379 2.728726 3.028567 1.000403 3685 4466
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 25.035088 3.432712 19.360955 32.834864 1.000438 5674 5556
##
## 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(ls2506_FLW_brm1))
qqline(resid(ls2506_FLW_brm1))
normal enough
Calculate R-squared
rss_fitted <- sum(residuals(ls2506_FLW_brm1)^2)
rss_constant <- sum((lsdf_2506_FL_W$W_G - mean(lsdf_2506_FL_W$W_G))^2)
approx_r_squared <- 1 - rss_fitted / rss_constant
approx_r_squared ## [1] 0.6607119smaller than the nls r-squared
Calculate the predicted W_G values based on this NLS model
ls2506_FLW_brm1_x <- c(16:34) #rounded actual min and max TL values
ls2506_FLW_brm1_predic <- data.frame(ls2506_FLW_brm1_x)
ls2506_FLW_brm1_f <- function(ls2506_FLW_brm1_x) {
#formula using calculated parameters
return(0.026165*ls2506_FLW_brm1_x^2.865512)
}
ls2506_FLW_brm1_predic$ls2506_FLW_brm1_y <- sapply(ls2506_FLW_brm1_predic$ls2506_FLW_brm1_x, ls2506_FLW_brm1_f)
View(ls2506_FLW_brm1_predic)Plot
ls2506_FLW_brm1_plottemp <- ls2506_FLW_nls_plot +
geom_smooth(aes(x = ls2506_FLW_brm1_x, y = ls2506_FLW_brm1_y),
data = ls2506_FLW_brm1_predic,
method = "loess",
se = F,
color = safe_colorblind_palette[4],
linewidth = 0.75,
alpha = 0.75)
# add label with equation
ls2506_FLW_brm1_plot <- ls2506_FLW_brm1_plottemp +
geom_text(x = 16, y = 470,
label = expression(W~'(g)' == 2.617%*%10^{-02} * TL~'(cm)'^{2.866}~"," ~ R^{2} == 0.656 ~ '(BRMS)'),
colour = safe_colorblind_palette[4],
hjust = 0)
ls2506_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.040314-0.014668## [1] 0.025646#B UCI-LCI
3.025578-2.728763## [1] 0.296815set priors and run model
#set priors
ls2506_FLW_prior2 <- prior(normal(0.026165, 0.025646) ,nlpar = "A")+prior(normal(2.865512, 0.296815),nlpar = "B")
ls2506_FLW_brm2 <- brm(bf(W_G~A*TL_CM^B, A + B ~ 1, nl = TRUE),data = lsdf_2506_FL_W,
prior = ls2506_FLW_prior2,
warmup = 1000,
iter = 5000,
chains = 4,
cores = 4,
control = list(adapt_delta = 0.8))
print(summary(ls2506_FLW_brm2), digits = 6)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: W_G ~ A * TL_CM^B
## A ~ 1
## B ~ 1
## Data: lsdf_2506_FL_W (Number of observations: 30)
## 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.018964 0.005931 0.010177 0.033242 1.000505 4003 4054
## B_Intercept 2.908418 0.086977 2.731779 3.075793 1.000488 4005 4049
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 19.780966 2.796125 15.251473 26.063945 1.000496 5110 5402
##
## 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(ls2506_FLW_brm2))
negative skew
Calculate R-squared
rss_fitted <- sum(residuals(ls2506_FLW_brm2)^2)
rss_constant <- sum((lsdf_2506_FL_W$W_G - mean(lsdf_2506_FL_W$W_G))^2)
approx_r_squared <- 1 - rss_fitted / rss_constant
approx_r_squared ## [1] 0.7868583better fit than prevoius brms but more skew of residuals
Calculate the predicted W_G values based on this NLS model
ls2506_FLW_brm2_x <- c(16:34) #rounded actual min and max TL values
ls2506_FLW_brm2_predic <- data.frame(ls2506_FLW_brm2_x)
ls2506_FLW_brm2_f <- function(ls2506_FLW_brm2_x) {
#formula using calculated parameters
return(0.018463*ls2506_FLW_brm2_x^2.917669)
}
ls2506_FLW_brm2_predic$ls2506_FLW_brm2_y <- sapply(ls2506_FLW_brm2_predic$ls2506_FLW_brm2_x, ls2506_FLW_brm2_f)
View(ls2506_FLW_brm2_predic)Plot
ls2506_FLW_brm2_plottemp <- ls2506_FLW_brm1_plot +
geom_smooth(aes(x = ls2506_FLW_brm2_x, y = ls2506_FLW_brm2_y),
data = ls2506_FLW_brm2_predic,
method = "loess",
se = F,
color = safe_colorblind_palette[5],
linewidth = 0.75)
# add label with equation
ls2506_FLW_brm2_plot <- ls2506_FLW_brm2_plottemp +
geom_text(x = 16, y = 430,
label = expression(W~'(g)' == 1.846%*%10^{-02} * TL~'(cm)'^{2.918}~"," ~ R^{2} == 0.788 ~ '(BRMS)'),
colour = safe_colorblind_palette[5],
hjust = 0)
ls2506_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
Plot the model
ls2506_FLW_nls_plot_title <- expression(italic("L. synagris")~"FL-W Relationship (n = 30)")
ls2506_FLW_nls_plottemp <- ggplot(lsdf_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 = ls2506_FLW_nls_plot_title) +
theme_classic() +
geom_smooth(aes(x = ls2506_FLW_nls_x, y = ls2506_FLW_nls_y),
data = ls2506_FLW_nls_predic,
method = "loess",
se = F,
color = "#A2FC3CFF",
linewidth = 0.75,
alpha = 0.75) +
ylim(0, 1600) +
xlim(12, 48)
# add label with equation
ls2506_FLW_nls_plot <- ls2506_FLW_nls_plottemp +
geom_label(x = 12, y = 1500,
label = expression(W~'(g)' == 2.487%*%10^{-02} * FL~'(cm)'^{2.871}~"," ~ R^{2} == 0.975 ~ '(NLS)'),
fill = "#c6ef97",
hjust = 0)
ls2506_FLW_nls_plot
Save plot
ggsave("ls2506_FLW_nls_plot.png", plot = ls2506_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
lsdf_2506_TL_W <- read.csv("ls2506df.csv")
lsdf_2506_TL_W$catch_month <- factor(lsdf_2506_TL_W$catch_month)
lsdf_2506_TL_W$year <- factor(lsdf_2506_TL_W$year)
lsdf_2506_TL_W$species <- factor(lsdf_2506_TL_W$species)
lsdf_2506_TL_W$common_name <- factor(lsdf_2506_TL_W$common_name)
lsdf_2506_TL_W$frozen <- factor(lsdf_2506_TL_W$frozen)
lsdf_2506_TL_W$catch_method <- factor(lsdf_2506_TL_W$catch_method)
lsdf_2506_TL_W$bait <- factor(lsdf_2506_TL_W$bait)
lsdf_2506_TL_W$area <- factor(lsdf_2506_TL_W$area)
lsdf_2506_TL_W$gutted <- factor(lsdf_2506_TL_W$gutted)
lsdf_2506_TL_W$sex <- factor(lsdf_2506_TL_W$sex)
lsdf_2506_TL_W$maturity <- factor(lsdf_2506_TL_W$maturity)
# Remove gutted samples
lsdf_2506_TL_W <- subset(filter(lsdf_2506_TL_W, gutted == "N"))
# remove NA rows
lsdf_2506_TL_W <- lsdf_2506_TL_W%>%
drop_na(TL_CM)
View(lsdf_2506_TL_W)
summary(lsdf_2506_TL_W)## X catch_date catch_month landing_time
## Min. : 1.00 Length:30 Apr: 1 Length:30
## 1st Qu.: 8.25 Class :character Aug: 3 Class :character
## Median :15.50 Mode :character Jul: 7 Mode :character
## Mean :15.50 Jun: 4
## 3rd Qu.:22.75 May:15
## Max. :30.00
##
## dissection_date dissection_time year month
## Length:30 Length:30 2023:11 Length:30
## Class :character Class :character 2024:19 Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
## observer species common_name fish_code
## Length:30 Lutjanus synagris:30 Lane Snapper:30 Length:30
## Class :character Class :character
## Mode :character Mode :character
##
##
##
##
## frozen catch_method bait area SL_CM
## N:19 :11 :11 :11 Min. :13.70
## Y:11 Lobster trap:19 Cowskin:19 B3 : 1 1st Qu.:19.00
## B3/4 : 2 Median :23.40
## B4 : 9 Mean :22.55
## B5 : 2 3rd Qu.:26.60
## D3 : 2 Max. :28.50
## LEAYLE: 3 NA's :11
## FL_CM TL_CM HH_CM ED_CM
## Min. :16.10 Min. :16.90 Min. :2.900 Min. :0.800
## 1st Qu.:23.82 1st Qu.:25.60 1st Qu.:4.200 1st Qu.:1.100
## Median :26.65 Median :28.75 Median :5.500 Median :1.250
## Mean :26.28 Mean :28.05 Mean :5.292 Mean :1.283
## 3rd Qu.:30.35 3rd Qu.:32.40 3rd Qu.:6.175 3rd Qu.:1.550
## Max. :33.50 Max. :35.10 Max. :7.000 Max. :1.900
## NA's :12 NA's :12
## W_G finclip otoliths gutted
## Min. : 65.7 Length:30 Length:30 N:30
## 1st Qu.:212.2 Class :character Class :character
## Median :315.2 Mode :character Mode :character
## Mean :319.8
## 3rd Qu.:454.9
## Max. :563.2
##
## gonads_present sex maturity G1L_CM G2L_CM
## Length:30 F:20 :11 Min. :2.000 Min. :1.200
## Class :character M:10 A: 1 1st Qu.:3.825 1st Qu.:4.000
## Mode :character B: 3 Median :4.750 Median :4.900
## C:15 Mean :4.630 Mean :4.910
## 3rd Qu.:5.375 3rd Qu.:5.975
## Max. :7.700 Max. :8.200
##
## GW_G stomach._sample
## Length:30 Length:30
## Class :character Class :character
## Mode :character Mode :character
##
##
##
## First look at the data with a scatterplot
fill = sex
ggplot(lsdf_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 synagris, n = 30") +
scale_fill_viridis(option = "D", discrete = TRUE)
Linear model
ls2506_TLW_lm1 <- lm(TL_CM ~ W_G, data = lsdf_2506_TL_W)
summary(ls2506_TLW_lm1)##
## Call:
## lm(formula = TL_CM ~ W_G, data = lsdf_2506_TL_W)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.6855 -1.0473 0.3523 0.8140 1.4570
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 17.398126 0.493001 35.29 <2e-16 ***
## W_G 0.033293 0.001397 23.84 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.142 on 28 degrees of freedom
## Multiple R-squared: 0.953, Adjusted R-squared: 0.9513
## F-statistic: 568.1 on 1 and 28 DF, p-value: < 2.2e-16Scatterplot the residuals against the fitted values to check homogeneity of variance
plot(resid(ls2506_TLW_lm1) ~ fitted(ls2506_TLW_lm1))
looks nonlinear
Plot a histogram to check normality of residuals
qqnorm(resid(ls2506_TLW_lm1))
qqline(resid(ls2506_TLW_lm1))
not at all normal
Plot the model
ls2506_TLW_lm1_plottemp <- ggplot(lsdf_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 synagris, n = 30") +
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
ls2506_TLW_lm1_plot <- ls2506_TLW_lm1_plottemp +
geom_text(x = 17, y = 550,
label = expression(W~'(g)' == 0.033 * TL~'(cm)' + 17.398~"," ~ R^{2} == 0.951 ~ "(Linear)"),
colour = safe_colorblind_palette[1],
hjust = 0)
ls2506_TLW_lm1_plot
NLS
ls2506_TLW_nls <- nls_multstart(W_G~A*TL_CM^B,
data = lsdf_2506_TL_W,
iter = 500,
start_lower = c(A=0, B=2),
start_upper = c(A=0.5, B=3.5))
print(summary(ls2506_TLW_nls), digits=5)##
## Formula: W_G ~ A * TL_CM^B
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## A 0.0158534 0.0048603 3.2618 0.00291 **
## B 2.9472123 0.0888054 33.1873 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 18.878 on 28 degrees of freedom
##
## Number of iterations to convergence: 27
## Achieved convergence tolerance: 1.4901e-08Plot a histogram of the residuals to check normality of residuals
qqnorm(resid(ls2506_TLW_nls))
qqline(resid(ls2506_TLW_nls))
fine
Calculate R-squared
rss_fitted <- sum(residuals(ls2506_TLW_nls)^2)
rss_constant <- sum((lsdf_2506_TL_W$W_G - mean(lsdf_2506_TL_W$W_G))^2)
approx_r_squared <- 1 - rss_fitted / rss_constant
approx_r_squared## [1] 0.9850624Calculate AIC
AIC(ls2506_TLW_nls)## [1] 265.3469Calculate the predicted W_kg values based on this NLS model
ls2506_TLW_nls_x <- c(16:36) #rounded actual min and max TL_MM values
ls2506_TLW_nls_predic <- data.frame(ls2506_TLW_nls_x)
ls2506_TLW_nls_f <- function(ls2506_TLW_nls_x) {
#formula using calculated parameters
return(0.01585*ls2506_TLW_nls_x^2.94721)
}
ls2506_TLW_nls_predic$ls2506_TLW_nls_y <- sapply(ls2506_TLW_nls_predic$ls2506_TLW_nls_x, ls2506_TLW_nls_f)
View(ls2506_TLW_nls_predic)Plot the model on the same graph as the other model
ls2506_TLW_nls_plottemp <- ls2506_TLW_lm1_plot +
geom_smooth(aes(x = ls2506_TLW_nls_x, y = ls2506_TLW_nls_y),
data = ls2506_TLW_nls_predic,
method = "loess",
se = F,
color = safe_colorblind_palette[2],
linewidth = 0.75)
# add label with equation
ls2506_TLW_nls_plot <- ls2506_TLW_nls_plottemp +
geom_text(x = 17, y = 500,
label = expression(W~'(g)' == 1.585%*%10^{-02} * TL~'(cm)'^{2.947}~"," ~ R^{2} == 0.985 ~ '(NLS)'),
colour = safe_colorblind_palette[2],
hjust = 0)
ls2506_TLW_nls_plot
Brms: with nls parameters to start
set priors and run model
#set priors
ls2506_TLW_prior1 <- prior(normal(0.01585, 0.005) ,nlpar = "A")+prior(normal(2.94721, 1),nlpar = "B")
ls2506_TLW_brm1 <- brm(bf(W_G~A*TL_CM^B, A + B ~ 1, nl = TRUE),data = lsdf_2506_TL_W,
prior = ls2506_TLW_prior1,
warmup = 1000,
iter = 5000,
chains = 4,
cores = 4,
control = list(adapt_delta = 0.8))
print(summary(ls2506_TLW_brm1), digits = 6)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: W_G ~ A * TL_CM^B
## A ~ 1
## B ~ 1
## Data: lsdf_2506_TL_W (Number of observations: 30)
## 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.016381 0.003477 0.010014 0.023705 1.000512 4502 4828
## B_Intercept 2.944423 0.063202 2.830664 3.079700 1.000455 4488 4797
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 19.517001 2.725339 15.077434 25.627983 1.000245 5220 5545
##
## 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(ls2506_TLW_brm1))
negative skew
Calculate R-squared
rss_fitted <- sum(residuals(ls2506_TLW_brm1)^2)
rss_constant <- sum((lsdf_2506_TL_W$W_G - mean(lsdf_2506_TL_W$W_G))^2)
approx_r_squared <- 1 - rss_fitted / rss_constant
approx_r_squared ## [1] 0.7950874Calculate the predicted W_G values based on this NLS model
ls2506_TLW_brm1_x <- c(16:36) #rounded actual min and max TL values
ls2506_TLW_brm1_predic <- data.frame(ls2506_TLW_brm1_x)
ls2506_TLW_brm1_f <- function(ls2506_TLW_brm1_x) {
#formula using calculated parameters
return(0.016379*ls2506_TLW_brm1_x^2.944415)
}
ls2506_TLW_brm1_predic$ls2506_TLW_brm1_y <- sapply(ls2506_TLW_brm1_predic$ls2506_TLW_brm1_x, ls2506_TLW_brm1_f)
View(ls2506_TLW_brm1_predic)Plot
ls2506_TLW_brm1_plottemp <- ls2506_TLW_nls_plot +
geom_smooth(aes(x = ls2506_TLW_brm1_x, y = ls2506_TLW_brm1_y),
data = ls2506_TLW_brm1_predic,
method = "loess",
se = F,
color = safe_colorblind_palette[4],
linewidth = 0.75)
# add label with equation
ls2506_TLW_brm1_plot <- ls2506_TLW_brm1_plottemp +
geom_text(x = 17, y = 450,
label = expression(W~'(g)' == 1.638%*%10^{-02} * TL~'(cm)'^{2.944}~"," ~ R^{2} == 0.794 ~ '(BRMS)'),
colour = safe_colorblind_palette[4],
hjust = 0)
ls2506_TLW_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.023590-0.010017## [1] 0.013573#B UCI-LCI
3.079806-2.831170 ## [1] 0.248636set priors and run model
#set priors
ls2506_TLW_prior2 <- prior(normal(0.016379, 0.013573) ,nlpar = "A")+prior(normal(2.944415, 0.248636),nlpar = "B")
ls2506_TLW_brm2 <- brm(bf(W_G~A*TL_CM^B, A + B ~ 1, nl = TRUE), data = lsdf_2506_TL_W,
prior = ls2506_TLW_prior2,
warmup = 1000,
iter = 5000,
chains = 4,
cores = 4,
control = list(adapt_delta = 0.8))
print(summary(ls2506_TLW_brm2), digits = 6)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: W_G ~ A * TL_CM^B
## A ~ 1
## B ~ 1
## Data: lsdf_2506_TL_W (Number of observations: 30)
## 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.017407 0.004825 0.009833 0.028302 1.003166 1425 2211
## B_Intercept 2.931010 0.080180 2.778752 3.085367 1.003161 1421 2505
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 19.624774 2.897507 14.305362 25.874157 1.011788 438 97
##
## 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(ls2506_TLW_brm2))
negative skew
Calculate R-squared
rss_fitted <- sum(residuals(ls2506_TLW_brm2)^2)
rss_constant <- sum((lsdf_2506_TL_W$W_G - mean(lsdf_2506_TL_W$W_G))^2)
approx_r_squared <- 1 - rss_fitted / rss_constant
approx_r_squared ## [1] 0.7898475Calculate the predicted W_G values based on this NLS model
ls2506_TLW_brm2_x <- c(16:36) #rounded actual min and max TL values
ls2506_TLW_brm2_predic <- data.frame(ls2506_TLW_brm2_x)
ls2506_TLW_brm2_f <- function(ls2506_TLW_brm2_x) {
#formula using calculated parameters
return(0.017708*ls2506_TLW_brm2_x^2.926179)
}
ls2506_TLW_brm2_predic$ls2506_TLW_brm2_y <- sapply(ls2506_TLW_brm2_predic$ls2506_TLW_brm2_x, ls2506_TLW_brm2_f)
View(ls2506_TLW_brm2_predic)Plot
ls2506_TLW_brm2_plottemp <- ls2506_TLW_brm1_plot +
geom_smooth(aes(x = ls2506_TLW_brm2_x, y = ls2506_TLW_brm2_y),
data = ls2506_TLW_brm2_predic,
method = "loess",
se = F,
color = safe_colorblind_palette[5],
linewidth = 0.75,
linetype = 2)
# add label with equation
ls2506_TLW_brm2_plot <- ls2506_TLW_brm2_plottemp +
geom_text(x = 16, y = 400,
label = expression(W~'(g)' == 1.771%*%10^{-02} * TL~'(cm)'^{2.926}~"," ~ R^{2} == "0.790" ~ '(BRMS)'),
colour = safe_colorblind_palette[5],
hjust = 0)
ls2506_TLW_brm2_plot
since the assumption of normality was not really met by these models, I will do a data transformation
Data transformation
# Find optimal lambda for Box-Cox transformation
bc_TLW <- boxcox(TL_CM ~ W_G, data = lsdf_2506_TL_W)
optimal_lambda <- bc_TLW$x[which.max(bc_TLW$y)]
optimal_lambda## [1] 2#lambda = 2 is a squared transformation
TL_CMsquared <- (lsdf_2506_TL_W$TL_CM)^2
lsdf_2506_TL_W$TL_CMsquared <- TL_CMsquared
W_Gsquared <- (lsdf_2506_TL_W$W_G)^2
lsdf_2506_TL_W$W_Gsquared <- W_Gsquared
summary(lsdf_2506_TL_W)## X catch_date catch_month landing_time
## Min. : 1.00 Length:30 Apr: 1 Length:30
## 1st Qu.: 8.25 Class :character Aug: 3 Class :character
## Median :15.50 Mode :character Jul: 7 Mode :character
## Mean :15.50 Jun: 4
## 3rd Qu.:22.75 May:15
## Max. :30.00
##
## dissection_date dissection_time year month
## Length:30 Length:30 2023:11 Length:30
## Class :character Class :character 2024:19 Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
## observer species common_name fish_code
## Length:30 Lutjanus synagris:30 Lane Snapper:30 Length:30
## Class :character Class :character
## Mode :character Mode :character
##
##
##
##
## frozen catch_method bait area SL_CM
## N:19 :11 :11 :11 Min. :13.70
## Y:11 Lobster trap:19 Cowskin:19 B3 : 1 1st Qu.:19.00
## B3/4 : 2 Median :23.40
## B4 : 9 Mean :22.55
## B5 : 2 3rd Qu.:26.60
## D3 : 2 Max. :28.50
## LEAYLE: 3 NA's :11
## FL_CM TL_CM HH_CM ED_CM
## Min. :16.10 Min. :16.90 Min. :2.900 Min. :0.800
## 1st Qu.:23.82 1st Qu.:25.60 1st Qu.:4.200 1st Qu.:1.100
## Median :26.65 Median :28.75 Median :5.500 Median :1.250
## Mean :26.28 Mean :28.05 Mean :5.292 Mean :1.283
## 3rd Qu.:30.35 3rd Qu.:32.40 3rd Qu.:6.175 3rd Qu.:1.550
## Max. :33.50 Max. :35.10 Max. :7.000 Max. :1.900
## NA's :12 NA's :12
## W_G finclip otoliths gutted
## Min. : 65.7 Length:30 Length:30 N:30
## 1st Qu.:212.2 Class :character Class :character
## Median :315.2 Mode :character Mode :character
## Mean :319.8
## 3rd Qu.:454.9
## Max. :563.2
##
## gonads_present sex maturity G1L_CM G2L_CM
## Length:30 F:20 :11 Min. :2.000 Min. :1.200
## Class :character M:10 A: 1 1st Qu.:3.825 1st Qu.:4.000
## Mode :character B: 3 Median :4.750 Median :4.900
## C:15 Mean :4.630 Mean :4.910
## 3rd Qu.:5.375 3rd Qu.:5.975
## Max. :7.700 Max. :8.200
##
## GW_G stomach._sample TL_CMsquared W_Gsquared
## Length:30 Length:30 Min. : 285.6 Min. : 4316
## Class :character Class :character 1st Qu.: 655.4 1st Qu.: 45082
## Mode :character Mode :character Median : 826.6 Median : 99493
## Mean : 812.5 Mean :124565
## 3rd Qu.:1049.8 3rd Qu.:207157
## Max. :1232.0 Max. :317194
## Linear model
ls2506_TLW_lmsquare <- lm(TL_CMsquared ~ W_Gsquared, data = lsdf_2506_TL_W)
summary(ls2506_TLW_lmsquare)##
## Call:
## lm(formula = TL_CMsquared ~ W_Gsquared, data = lsdf_2506_TL_W)
##
## Residuals:
## Min 1Q Median 3Q Max
## -211.82 -57.62 43.81 66.04 108.46
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.861e+02 2.741e+01 17.73 < 2e-16 ***
## W_Gsquared 2.620e-03 1.730e-04 15.15 5.09e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 92.83 on 28 degrees of freedom
## Multiple R-squared: 0.8913, Adjusted R-squared: 0.8874
## F-statistic: 229.5 on 1 and 28 DF, p-value: 5.09e-15Scatterplot the residuals against the fitted values to check homogeneity of variance
plot(resid(ls2506_TLW_lmsquare) ~ fitted(ls2506_TLW_lmsquare))
looks nonlinear
Plot a histogram to check normality of residuals
qqnorm(resid(ls2506_TLW_lmsquare))
qqline(resid(ls2506_TLW_lmsquare))
more normal than before
Plot the model
ls2506_TLW_lmsquare_plottemp <- ggplot(lsdf_2506_TL_W, aes(x = TL_CMsquared, y = W_Gsquared)) +
geom_point(aes(TL_CMsquared, W_Gsquared),
color = "#242124",
alpha = 0.8) +
labs(x = expression("(Total Length in cm)"^{2}),
y = expression("(Weight in grams)"^{2}),
title = "Lutjanis synagris, n = 30") +
theme_classic() +
geom_smooth(method = lm,
aes(TL_CMsquared, W_Gsquared),
linewidth = 0.5,
se = FALSE,
color = safe_colorblind_palette[1])
# add label with equation
ls2506_TLW_lmsquare_plot <- ls2506_TLW_lmsquare_plottemp +
geom_text(x = 300, y = 300000,
label = expression(W~'(g)' == 0.033 * TL~'(cm)' + 17.398~"," ~ R^{2} == 0.951 ~ "(Linear)"),
colour = safe_colorblind_palette[1],
hjust = 0)
ls2506_TLW_lmsquare_plot
NLS
ls2506_TLW_nlssquare <- nls_multstart(W_Gsquared~A*TL_CMsquared^B,
data = lsdf_2506_TL_W,
iter = 500,
start_lower = c(A=0, B=2),
start_upper = c(A=0.5, B=3.5))
summary(ls2506_TLW_nlssquare)##
## Formula: W_Gsquared ~ A * TL_CMsquared^B
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## A 0.0004398 0.0004051 1.086 0.287
## B 2.8671074 0.1317082 21.769 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 16650 on 28 degrees of freedom
##
## Number of iterations to convergence: 42
## Achieved convergence tolerance: 1.49e-08Plot a histogram of the residuals to check normality of residuals
qqnorm(resid(ls2506_TLW_nlssquare))
qqline(resid(ls2506_TLW_nlssquare))
heavy tails but okay
Calculate R-squared
rss_fitted <- sum(residuals(ls2506_TLW_nlssquare)^2)
rss_constant <- sum((lsdf_2506_TL_W$W_Gsquared- mean(lsdf_2506_TL_W$W_Gsquared))^2)
approx_r_squared <- 1 - rss_fitted / rss_constant
approx_r_squared## [1] 0.9730526Calculate AIC
AIC(ls2506_TLW_nlssquare)## [1] 672.2808Calculate the predicted W_kg values based on this NLS model
ls2506_TLW_nlssquare_x <- c(285:1232) #rounded actual min and max TL_MM values
ls2506_TLW_nlssquare_predic <- data.frame(ls2506_TLW_nlssquare_x)
ls2506_TLW_nlssquare_f <- function(ls2506_TLW_nlssquare_x) {
#formula using calculated parameters
return(0.0004398*ls2506_TLW_nlssquare_x^2.8671073)
}
ls2506_TLW_nlssquare_predic$ls2506_TLW_nlssquare_y <- sapply(ls2506_TLW_nlssquare_predic$ls2506_TLW_nlssquare_x, ls2506_TLW_nlssquare_f)
View(ls2506_TLW_nlssquare_predic)Plot the model on the same graph as the other model
ls2506_TLW_nlssquare_plottemp <- ls2506_TLW_lmsquare_plot +
geom_smooth(aes(x = ls2506_TLW_nlssquare_x, y = ls2506_TLW_nlssquare_y),
data = ls2506_TLW_nlssquare_predic,
method = "loess",
se = F,
color = safe_colorblind_palette[2],
linewidth = 0.75)
# add label with equation
ls2506_TLW_nlssquare_plot <- ls2506_TLW_nlssquare_plottemp +
geom_text(x = 300, y = 270000,
label = expression(W^{2}~'(g)' == 4.398%*%10^{-04} * TL^{2}~'(cm)'^{2.867}~"," ~ R^{2} == 0.973 ~ '(NLS)'),
colour = safe_colorblind_palette[2],
hjust = 0)
ls2506_TLW_nlssquare_plot
##Final decision: NLS with variables transformed by squaring
Plot the model
ls2506_sqTLW_nls_plot_title <- expression(italic("L. synagris") ~ "TL-W Relationship (n = 30)")
ls2506_sqTLW_nls_plottemp <- ggplot(lsdf_2506_FL_W, aes(x = TL_CMsquared, y = W_Gsquared)) +
geom_point(aes(TL_CMsquared, W_Gsquared),
color = "#242124",
alpha = 0.8) +
labs(x = expression("(Total Length in cm)"^{2}),
y = expression("(Weight in grams)"^{2}), # Corrected closing parenthesis
title = ls2506_sqTLW_nls_plot_title) +
theme_classic() +
geom_smooth(aes(x = ls2506_TLW_nlssquare_x, y = ls2506_TLW_nlssquare_y),
data = ls2506_TLW_nlssquare_predic,
method = "loess",
se = FALSE,
color = "#61b303",
linewidth = 0.75)
# add label with equation
ls2506_sqTLW_nls_plot <- ls2506_sqTLW_nls_plottemp +
geom_text(x = 300, y = 300000,
label = expression(W^{2}~'(g)' == 4.398%*%10^{-04} * TL^{2}~'(cm)'^{2.867}~"," ~ R^{2} == 0.973 ~ '(NLS)'),
colour = "#61b303",
hjust = 0)
ls2506_sqTLW_nls_plot
Standard length (SL_CM) - Weight (W_G) relationship
Load data file
lsdf_2506_SL_W <- read.csv("ls2506df.csv")
lsdf_2506_SL_W$catch_month <- factor(lsdf_2506_SL_W$catch_month)
lsdf_2506_SL_W$year <- factor(lsdf_2506_SL_W$year)
lsdf_2506_SL_W$species <- factor(lsdf_2506_SL_W$species)
lsdf_2506_SL_W$common_name <- factor(lsdf_2506_SL_W$common_name)
lsdf_2506_SL_W$frozen <- factor(lsdf_2506_SL_W$frozen)
lsdf_2506_SL_W$catch_method <- factor(lsdf_2506_SL_W$catch_method)
lsdf_2506_SL_W$bait <- factor(lsdf_2506_SL_W$bait)
lsdf_2506_SL_W$area <- factor(lsdf_2506_SL_W$area)
lsdf_2506_SL_W$gutted <- factor(lsdf_2506_SL_W$gutted)
lsdf_2506_SL_W$sex <- factor(lsdf_2506_SL_W$sex)
lsdf_2506_SL_W$maturity <- factor(lsdf_2506_SL_W$maturity)
# Remove gutted samples
lsdf_2506_SL_W <- subset(filter(lsdf_2506_SL_W, gutted == "N"))
# remove NA rows
lsdf_2506_SL_W <- lsdf_2506_SL_W%>%
drop_na(SL_CM)
View(lsdf_2506_SL_W)
summary(lsdf_2506_SL_W)## X catch_date catch_month landing_time
## Min. : 1.0 Length:19 Apr: 1 Length:19
## 1st Qu.: 5.5 Class :character Aug: 0 Class :character
## Median :10.0 Mode :character Jul: 0 Mode :character
## Mean :10.0 Jun: 3
## 3rd Qu.:14.5 May:15
## Max. :19.0
##
## dissection_date dissection_time year month
## Length:19 Length:19 2023: 0 Length:19
## Class :character Class :character 2024:19 Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
## observer species common_name fish_code
## Length:19 Lutjanus synagris:19 Lane Snapper:19 Length:19
## Class :character Class :character
## Mode :character Mode :character
##
##
##
##
## frozen catch_method bait area SL_CM
## N:19 : 0 : 0 :0 Min. :13.70
## Y: 0 Lobster trap:19 Cowskin:19 B3 :1 1st Qu.:19.00
## B3/4 :2 Median :23.40
## B4 :9 Mean :22.55
## B5 :2 3rd Qu.:26.60
## D3 :2 Max. :28.50
## LEAYLE:3
## FL_CM TL_CM HH_CM ED_CM
## Min. :16.10 Min. :16.90 Min. :2.900 Min. :0.800
## 1st Qu.:22.35 1st Qu.:23.65 1st Qu.:4.200 1st Qu.:1.100
## Median :26.70 Median :28.70 Median :5.500 Median :1.250
## Mean :25.84 Mean :27.53 Mean :5.292 Mean :1.283
## 3rd Qu.:30.60 3rd Qu.:32.55 3rd Qu.:6.175 3rd Qu.:1.550
## Max. :32.20 Max. :34.50 Max. :7.000 Max. :1.900
## NA's :1 NA's :1
## W_G finclip otoliths gutted
## Min. : 65.7 Length:19 Length:19 N:19
## 1st Qu.:169.5 Class :character Class :character
## Median :303.3 Mode :character Mode :character
## Mean :309.0
## 3rd Qu.:472.3
## Max. :562.9
##
## gonads_present sex maturity G1L_CM G2L_CM
## Length:19 F:11 : 0 Min. :2.500 Min. :1.200
## Class :character M: 8 A: 1 1st Qu.:3.250 1st Qu.:3.600
## Mode :character B: 3 Median :4.700 Median :4.600
## C:15 Mean :4.368 Mean :4.468
## 3rd Qu.:5.250 3rd Qu.:5.550
## Max. :6.100 Max. :7.200
##
## GW_G stomach._sample
## Length:19 Length:19
## Class :character Class :character
## Mode :character Mode :character
##
##
##
## First look at the data with a scatterplot
fill = sex
ggplot(lsdf_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 synagris, n = 19") +
scale_fill_viridis(option = "D", discrete = TRUE)
LM Linear model FL-W
ls2506_FLW_lm <- lm(SL_CM ~ W_G, data = lsdf_2506_SL_W)
summary(ls2506_FLW_lm)##
## Call:
## lm(formula = SL_CM ~ W_G, data = lsdf_2506_SL_W)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.8116 -0.8874 0.2845 0.6677 1.8869
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.611457 0.527214 25.82 4.45e-15 ***
## W_G 0.028922 0.001514 19.10 6.32e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.06 on 17 degrees of freedom
## Multiple R-squared: 0.9555, Adjusted R-squared: 0.9529
## F-statistic: 365 on 1 and 17 DF, p-value: 6.321e-13Scatterplot the residuals against the fitted values to check homogeneity of variance
plot(resid(ls2506_FLW_lm) ~ fitted(ls2506_FLW_lm))
looks like it indicates a nonlinear relationship
Plot a histogram to check normality of residuals
qqnorm(resid(ls2506_FLW_lm))
qqline(resid(ls2506_FLW_lm))
Plot the model
ls2506_SLW_lm_plottemp <- ggplot(lsdf_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 synagris, n = 19") +
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
ls2506_SLW_lm_plot <- ls2506_SLW_lm_plottemp +
geom_text(x = 13, y = 550,
label = expression(W~'(g)' == 0.029 * TL~'(cm)' + 13.611~"," ~ R^{2} == 0.9529 ~ "(Linear)"),
colour = safe_colorblind_palette[1],
hjust = 0)
ls2506_SLW_lm_plot
NLS Models
NLS 1
ls2506_SLW_nls <- nls_multstart(W_G~A*SL_CM^B,
data = lsdf_2506_SL_W,
iter = 500,
start_lower = c(A=0, B=2),
start_upper = c(A=0.5, B=3.5))
summary(ls2506_SLW_nls)##
## Formula: W_G ~ A * SL_CM^B
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## A 0.03347 0.01601 2.09 0.052 .
## B 2.89512 0.14680 19.72 3.76e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 24.49 on 17 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
#check residuals
qqnorm(resid(ls2506_SLW_nls))
qqline(resid(ls2506_SLW_nls))
okay
Calculate R-squared
rss_fitted <- sum(residuals(ls2506_SLW_nls)^2)
rss_constant <- sum((lsdf_2506_SL_W$W_G - mean(lsdf_2506_SL_W$W_G))^2)
approx_r_squared <- 1 - rss_fitted / rss_constant
approx_r_squared## [1] 0.9792107Calculate AIC
AIC(ls2506_SLW_nls)## [1] 179.3413Calculate the predicted W_g values based on this NLS model
ls2506_SLW_nls_x <- c(13:29) #rounded actual min and max FL_MM values
ls2506_SLW_nls_predic <- data.frame(ls2506_SLW_nls_x)
ls2506_SLW_nls_f <- function(ls2506_SLW_nls_x) {
#formula using calculated parameters
return(0.03347*ls2506_SLW_nls_x^2.89512)
}
ls2506_SLW_nls_predic$ls2506_SLW_nls_y <- sapply(ls2506_SLW_nls_predic$ls2506_SLW_nls_x, ls2506_SLW_nls_f)
View(ls2506_SLW_nls_predic)Plot the model on the same graph as the other model
ls2506_SLW_nls_plottemp <- ls2506_SLW_lm_plot +
geom_smooth(aes(x = ls2506_SLW_nls_x, y = ls2506_SLW_nls_y),
data = ls2506_SLW_nls_predic,
method = "loess",
se = F,
color = safe_colorblind_palette[2],
linewidth = 0.75,
alpha = 0.75)
# add label with equation
ls2506_SLW_nls_plot <- ls2506_SLW_nls_plottemp +
geom_text(x = 13, y = 510,
label = expression(W~'(g)' == 3.347%*%10^{-02} * TL~'(cm)'^{2.895}~"," ~ R^{2} == 0.979 ~ '(NLS)'),
colour = safe_colorblind_palette[2],
hjust = 0)
ls2506_SLW_nls_plot
##Final decision: NLS
Plot the model
ls2506_SLW_nls_plot_title <- expression(italic("L. synagris") ~ "SL-W Relationship (n = 19)")
ls2506_sqSLW_nls_plottemp <- ggplot(lsdf_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 = ls2506_SLW_nls_plot_title) +
theme_classic() +
geom_smooth(aes(x = ls2506_SLW_nls_x, y = ls2506_SLW_nls_y),
data = ls2506_SLW_nls_predic,
method = "loess",
se = F,
color = "#a9fc4c",
linewidth = 0.75,
alpha = 0.75)
# add label with equation
ls2506_sqSLW_nls_plot <- ls2506_sqSLW_nls_plottemp +
geom_text(x = 13, y = 520,
label = expression(W~'(g)' == 3.347%*%10^{-02} * TL~'(cm)'^{2.895}~"," ~ R^{2} == 0.979 ~ '(NLS)'),
colour = "#a9fc4c",
hjust = 0)
ls2506_sqTLW_nls_plot
Maturity
Load data file
lsdf_2506_maturity <- read.csv("ls2506df.csv")
lsdf_2506_maturity$catch_month <- factor(lsdf_2506_maturity$catch_month)
lsdf_2506_maturity$year <- factor(lsdf_2506_maturity$year)
lsdf_2506_maturity$species <- factor(lsdf_2506_maturity$species)
lsdf_2506_maturity$common_name <- factor(lsdf_2506_maturity$common_name)
lsdf_2506_maturity$frozen <- factor(lsdf_2506_maturity$frozen)
lsdf_2506_maturity$catch_method <- factor(lsdf_2506_maturity$catch_method)
lsdf_2506_maturity$bait <- factor(lsdf_2506_maturity$bait)
lsdf_2506_maturity$area <- factor(lsdf_2506_maturity$area)
lsdf_2506_maturity$gutted <- factor(lsdf_2506_maturity$gutted)
lsdf_2506_maturity$sex <- factor(lsdf_2506_maturity$sex)
lsdf_2506_maturity$maturity <- factor(lsdf_2506_maturity$maturity)
#drop the columns with neither sex nor maturity measurements
lsdf_2506_maturity <- lsdf_2506_maturity %>%
filter(maturity != "" | sex != "")calculate sex ratio
ls_totalfemale <- lsdf_2506_maturity %>%
filter(sex == "F") %>%
count()
ls_totalmale <- lsdf_2506_maturity %>%
filter(sex == "M") %>%
count()
ls_sexratio <- ls_totalfemale/ls_totalmale
ls_sexratio #female number## n
## 1 2Descriptive statistics
ls_sex_summary <- lsdf_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)
)
ls_sex_summary## # A tibble: 2 × 8
## sex count mean sd median min max IQR
## <fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 F 20 25.4 4.32 25.5 17.7 31.9 5.03
## 2 M 10 28.0 5.36 28.9 16.1 33.5 6.1lsdf_2506_maturity <- lsdf_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))
ls_maturity_summary <- lsdf_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)
)
ls_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 1 16.1 NA 16.1 16.1 16.1 0
## 2 mature 18 26.4 4.99 27.0 17.7 32.2 6.68
## 3 <NA> 11 27.0 3.56 26.6 22.6 33.5 3.70