---
title: "sizeMat: An R Package to Estimate Size at Sexual Maturity"
author: "Josymar Torrejon-Magallanes"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{sizeMat: An R Package to Estimate Size at Sexual Maturity}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.align = "center",
  fig.width = 7,
  fig.height = 5,
  warning = FALSE,
  message = FALSE
)

library(sizeMat)
```

# Introduction

The `sizeMat` package estimates size at sexual maturity in organisms, especially fish and invertebrates.

It includes tools for estimating:

- size at morphometric maturity, based on relative growth and classification of individuals into juvenile and adult groups;
- size at gonadal maturity, based on maturity stages;
- maturity ogives using frequentist or Bayesian logistic regression;
- uncertainty around model parameters and size at 50% maturity, usually denoted as $L_{50}$.

The size at sexual maturity is defined here as the length at which a randomly chosen individual has a 50% probability of being mature.

# Installation

The stable version can be installed from CRAN:

```{r install-cran, eval = FALSE}
install.packages("sizeMat")
```

The development version can be installed from GitHub:

```{r install-github, eval = FALSE}
# install.packages("devtools")
devtools::install_github("ejosymart/sizeMat")
```

# 1. Size at morphometric maturity

Morphometric maturity is estimated from the relative growth relationship between two allometric variables.

In this example, we use the `crabdata` dataset. This dataset contains allometric measurements and biological information for crabs of the species *Chionoecetes tanneri*.

```{r load-crabdata}
data(crabdata)

head(crabdata)
names(crabdata)
```

## 1.1 Classification of juveniles and adults

The function `classify_mature()` classifies individuals into two groups:

- juveniles = 0
- adults = 1

The classification is based on a principal components analysis using two allometric variables, followed by hierarchical clustering and discriminant analysis.

The argument `varNames` requires two allometric variables. The first variable is treated as the independent variable and the second as the dependent variable. The argument `varSex` identifies the sex variable. If `selectSex = NULL`, all individuals are used.

```{r classify-data}
classify_data <- classify_mature(
  crabdata,
  varNames = c("carapace_width", "chela_height"),
  varSex = "sex_category",
  selectSex = NULL,
  method = "ld"
)

classify_data_males <- classify_mature(
  crabdata,
  varNames = c("carapace_width", "chela_height"),
  varSex = "sex_category",
  selectSex = "m",
  method = "ld"
)

print(classify_data)
```

## 1.2 Base graphics plot for classified data

By default, `plot.classify()` uses base R graphics.

```{r plot-classify-base, fig.width = 7, fig.height = 5}
plot(
  classify_data,
  xlab = "Carapace width (mm)",
  ylab = "Chela height (mm)"
)
```

The graphic can be customized using colors, point symbols, line types, and point size.

```{r plot-classify-base-custom, fig.width = 7, fig.height = 5}
plot(
  classify_data,
  xlab = "Carapace width (mm)",
  ylab = "Chela height (mm)",
  col = c(2, 3),
  pch = c(5, 6),
  lty_lines = c(1, 2),
  lwd_lines = c(1, 3),
  cex = c(1, 2),
  main = "Classification"
)
```

## 1.3 ggplot2-style plot for classified data

A `ggplot2` style can be obtained with `gg_style = TRUE`.

```{r plot-classify-gg, fig.width = 7, fig.height = 5}
plot(
  classify_data,
  xlab = "Carapace width (mm)",
  ylab = "Chela height (mm)",
  col = c("steelblue", "firebrick"),
  pch = c(16, 17),
  gg_style = TRUE
)
```

Because the function returns a `ggplot` object when `gg_style = TRUE`, the plot can be further modified.

```{r plot-classify-gg-modified, fig.width = 7, fig.height = 5}
p_classify <- plot(
  classify_data,
  xlab = "Carapace width (mm)",
  ylab = "Chela height (mm)",
  col = c("steelblue", "firebrick"),
  pch = c(16, 17),
  gg_style = TRUE
)

p_classify + ggplot2::theme_bw()
```

# 2. Estimation of size at morphometric maturity

After individuals are classified, size at morphometric maturity can be estimated using `morph_mature()`.

The model is based on a logistic function:

$$
P_{CS} = \frac{1}{1 + e^{-(\hat{\beta}_{0} + \hat{\beta}_{1}X)}}
$$

where $P_{CS}$ is the probability of being mature at length $X$.

The size at 50% maturity is calculated as:

$$
L_{50} = -\frac{\hat{\beta}_{0}}{\hat{\beta}_{1}}
$$

The argument `method` defines the model type:

- `method = "fq"` uses frequentist logistic regression;
- `method = "bayes"` uses Bayesian logistic regression.

For vignette purposes, we use a relatively small number of iterations so that the document compiles quickly.

```{r morph-mature}
set.seed(123)

my_morph_fq <- morph_mature(
  classify_data,
  method = "fq",
  niter = 200
)

print(my_morph_fq)

my_morph_bayes <- morph_mature(
  classify_data,
  method = "bayes",
  niter = 200
)

print(my_morph_bayes)
```

## 2.1 Base graphics plot for morphometric maturity

By default, `plot.morphMat()` uses base R graphics. It produces histograms for the model parameters and the maturity ogive.

```{r plot-morph-base, fig.width = 7, fig.height = 5}
oldpar <- par(no.readonly = TRUE)
par(mfrow = c(2, 2))

plot(
  my_morph_fq,
  xlab = "Carapace width (mm)",
  ylab = "Proportion mature",
  col = c("blue", "red")
)

par(oldpar)
```

To plot only the maturity ogive, use `onlyOgive = TRUE`.

```{r plot-morph-base-ogive, fig.width = 7, fig.height = 5}
plot(
  my_morph_fq,
  xlab = "Carapace width (mm)",
  ylab = "Proportion mature",
  col = c("blue", "red"),
  onlyOgive = TRUE
)
```

## 2.2 ggplot2-style plots for morphometric maturity

A `ggplot2` maturity ogive can be obtained using `gg_style = TRUE`.

```{r plot-morph-gg-ogive, fig.width = 7, fig.height = 5}
plot(
  my_morph_fq,
  xlab = "Carapace width (mm)",
  ylab = "Proportion mature",
  col = c("steelblue", "firebrick"),
  onlyOgive = TRUE,
  gg_style = TRUE
)
```

When `onlyOgive = FALSE` and `gg_style = TRUE`, the function returns a list of independent `ggplot` objects.

```{r plot-morph-gg-list}
p_morph <- plot(
  my_morph_fq,
  xlab = "Carapace width (mm)",
  ylab = "Proportion mature",
  col = c("steelblue", "firebrick"),
  gg_style = TRUE
)

names(p_morph)
```

The individual plots can then be displayed separately.

```{r plot-morph-gg-A, fig.width = 7, fig.height = 5}
p_morph$A
```

```{r plot-morph-gg-B, fig.width = 7, fig.height = 5}
p_morph$B
```

```{r plot-morph-gg-L50, fig.width = 7, fig.height = 5}
p_morph$L50
```

```{r plot-morph-gg-ogive-from-list, fig.width = 7, fig.height = 5}
p_morph$ogive
```

The returned plots can also be modified using standard `ggplot2` syntax.

```{r plot-morph-gg-custom, fig.width = 7, fig.height = 5}
p_morph$ogive +
  ggplot2::theme_bw()
```

# 3. Size at gonadal maturity

Gonadal maturity is estimated using a maturity stage variable and an allometric variable such as total length.

In this example, we use the `matFish` dataset.

```{r load-matfish}
data(matFish)

head(matFish)
```

The dataset contains:

- `total_length`: total length in cm;
- `stage_mat`: gonadal maturation stage.

In this example, stage `I` is considered immature, whereas stages `II`, `III`, and `IV` are considered mature.

## 3.1 Estimation of size at gonadal maturity

The function `gonad_mature()` requires:

- `varNames`: the length variable and the maturity stage variable;
- `immName`: the immature stage or stages;
- `matName`: the mature stage or stages;
- `method`: the logistic regression approach;
- `niter`: the number of bootstrap or Bayesian iterations.

```{r gonad-mature}
set.seed(123)

my_gonad_fq <- gonad_mature(
  matFish,
  varNames = c("total_length", "stage_mat"),
  immName = "I",
  matName = c("II", "III", "IV"),
  method = "fq",
  niter = 200
)

print(my_gonad_fq)

my_gonad_bayes <- gonad_mature(
  matFish,
  varNames = c("total_length", "stage_mat"),
  immName = "I",
  matName = c("II", "III", "IV"),
  method = "bayes",
  niter = 200
)

print(my_gonad_bayes)
```

## 3.2 Base graphics plot for gonadal maturity

By default, `plot.gonadMat()` uses base R graphics.

```{r plot-gonad-base, fig.width = 7, fig.height = 5}
oldpar <- par(no.readonly = TRUE)
par(mfrow = c(2, 2))

plot(
  my_gonad_fq,
  xlab = "Total length (cm)",
  ylab = "Proportion mature",
  col = c("blue", "red")
)

par(oldpar)
```

To plot only the maturity ogive, use `onlyOgive = TRUE`.

```{r plot-gonad-base-ogive, fig.width = 7, fig.height = 5}
plot(
  my_gonad_fq,
  xlab = "Total length (cm)",
  ylab = "Proportion mature",
  col = c("blue", "red"),
  onlyOgive = TRUE
)
```

The legend position can be modified with `legendPosition`.

```{r plot-gonad-base-legend-position, fig.width = 7, fig.height = 5}
plot(
  my_gonad_fq,
  xlab = "Total length (cm)",
  ylab = "Proportion mature",
  col = c("blue", "red"),
  onlyOgive = TRUE,
  legendPosition = "bottomright"
)
```

The legend can also be removed.

```{r plot-gonad-base-no-legend, fig.width = 7, fig.height = 5}
plot(
  my_gonad_fq,
  xlab = "Total length (cm)",
  ylab = "Proportion mature",
  col = c("blue", "red"),
  onlyOgive = TRUE,
  showLegend = FALSE
)
```

## 3.3 ggplot2-style plots for gonadal maturity

A `ggplot2` maturity ogive can be obtained using `gg_style = TRUE`.

```{r plot-gonad-gg-ogive, fig.width = 7, fig.height = 5}
plot(
  my_gonad_fq,
  xlab = "Total length (cm)",
  ylab = "Proportion mature",
  col = c("steelblue", "firebrick"),
  onlyOgive = TRUE,
  gg_style = TRUE
)
```

The position of the $L_{50}$ and $R^2$ labels can also be changed.

```{r plot-gonad-gg-legend-position, fig.width = 7, fig.height = 5}
plot(
  my_gonad_fq,
  xlab = "Total length (cm)",
  ylab = "Proportion mature",
  col = c("steelblue", "firebrick"),
  onlyOgive = TRUE,
  gg_style = TRUE,
  legendPosition = "bottomright"
)
```

When `onlyOgive = FALSE` and `gg_style = TRUE`, the function returns a list of independent `ggplot` objects.

```{r plot-gonad-gg-list}
p_gonad <- plot(
  my_gonad_fq,
  xlab = "Total length (cm)",
  ylab = "Proportion mature",
  col = c("steelblue", "firebrick"),
  gg_style = TRUE
)

names(p_gonad)
```

The individual plots can be displayed separately.

```{r plot-gonad-gg-A, fig.width = 7, fig.height = 5}
p_gonad$A
```

```{r plot-gonad-gg-B, fig.width = 7, fig.height = 5}
p_gonad$B
```

```{r plot-gonad-gg-L50, fig.width = 7, fig.height = 5}
p_gonad$L50
```

```{r plot-gonad-gg-ogive-from-list, fig.width = 7, fig.height = 5}
p_gonad$ogive
```

The returned plots can be customized using `ggplot2`.

```{r plot-gonad-gg-custom, fig.width = 7, fig.height = 5}
p_gonad$ogive +
  ggplot2::theme_bw()
```

# References

Agostinho, C. S. (2000). Use of otoliths to estimate size at sexual maturity in fish. *Brazilian Archives of Biology and Technology*, 43(4).

Corgos, A. and Freire, J. (2006). Morphometric and gonad maturity in the spider crab *Maja brachydactyla*: a comparison of methods for estimating size at maturity in species with determinate growth. *ICES Journal of Marine Science*, 63(5), 851-859.

Roa, R., Ernst, B. and Tapia, F. (1999). Estimation of size at sexual maturity: an evaluation of analytical and resampling procedures. *Fishery Bulletin*, 97(3), 570-580.

Somerton, D. A. (1980). A computer technique for estimating the size of sexual maturity in crabs. *Canadian Journal of Fisheries and Aquatic Sciences*, 37(10), 1488-1494.