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:

install.packages("sizeMat")

The development version can be installed from GitHub:

# 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.

data(crabdata)

head(crabdata)
#>   year month carapace_width carapace_length chela_height chela_width
#> 1 1974     1            106             107         14.0          22
#> 2 1974     1            129             129         27.0          44
#> 3 1974     1            119             122         14.6          23
#> 4 1974     1            115             118         18.6          29
#> 5 1974     1             97              97         11.0          17
#> 6 1974     1             94              96         10.0          15
#>   sex_category
#> 1            m
#> 2            m
#> 3            m
#> 4            m
#> 5            m
#> 6            m
names(crabdata)
#> [1] "year"            "month"           "carapace_width"  "carapace_length"
#> [5] "chela_height"    "chela_width"     "sex_category"

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.

classify_data <- classify_mature(
  crabdata,
  varNames = c("carapace_width", "chela_height"),
  varSex = "sex_category",
  selectSex = NULL,
  method = "ld"
)
#> all individuals were used in the analysis

classify_data_males <- classify_mature(
  crabdata,
  varNames = c("carapace_width", "chela_height"),
  varSex = "sex_category",
  selectSex = "m",
  method = "ld"
)
#> only m-sex were used in the analysis

print(classify_data)
#> Number in juvenile group = 83 
#> 
#> Number in adult group = 140 
#> 
#> -------------------------------------------------------- 
#> 1) Linear regression for juveniles 
#> 
#> Call:
#> glm(formula = y ~ x, data = juv)
#> 
#> Coefficients:
#>              Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) -3.794687   0.497056  -7.634 3.93e-11 ***
#> x            0.161327   0.004701  34.314  < 2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for gaussian family taken to be 0.7320842)
#> 
#>     Null deviance: 921.306  on 82  degrees of freedom
#> Residual deviance:  59.299  on 81  degrees of freedom
#> AIC: 213.63
#> 
#> Number of Fisher Scoring iterations: 2
#> 
#> -------------------------------------------------------- 
#> 2) Linear regression for adults 
#> 
#> Call:
#> glm(formula = y ~ x, data = adt)
#> 
#> Coefficients:
#>               Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) -11.246726   1.199496  -9.376   <2e-16 ***
#> x             0.273837   0.008648  31.663   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for gaussian family taken to be 2.265729)
#> 
#>     Null deviance: 2584.24  on 139  degrees of freedom
#> Residual deviance:  312.67  on 138  degrees of freedom
#> AIC: 515.79
#> 
#> Number of Fisher Scoring iterations: 2
#> 
#> -------------------------------------------------------- 
#> 3) Difference between slopes (ANCOVA) 
#>               Estimate  Std. Error   t value     Pr(>|t|)
#> (Intercept) -3.7946869 0.757105677 -5.012097 1.109526e-06
#> x            0.1613275 0.007161179 22.528064 6.035478e-59
#> mature      -7.4520389 1.285219562 -5.798261 2.320729e-08
#> x:mature     0.1125093 0.010361046 10.858878 2.956242e-22
#> [1] "slopes are different"

1.2 Base graphics plot for classified data

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

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.

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.

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.

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.

set.seed(123)

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

print(my_morph_fq)
#> formula: Y = 1/1+exp-(A + B*X)
#>     Original Bootstrap (Median)
#> A   -20.753  -21.086           
#> B   0.1748   0.1762            
#> L50 118.7237 118.6442          
#> R2  0.7111   -

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

print(my_morph_bayes)
#> formula: Y = 1/1+exp-(A + B*X)
#>     Bootstrap (Median)
#> A             -20.3108
#> B               0.1724
#> L50           118.3878
#> R2              0.7111

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.

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")
)

#> Size at morphometric maturity = 118.6 
#> Confidence intervals = 116 - 121 
#> Rsquare = 0.71

par(oldpar)

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

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

#> Size at morphometric maturity = 118.6 
#> Confidence intervals = 116 - 121 
#> Rsquare = 0.71

2.2 ggplot2-style plots for morphometric maturity

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

plot(
  my_morph_fq,
  xlab = "Carapace width (mm)",
  ylab = "Proportion mature",
  col = c("steelblue", "firebrick"),
  onlyOgive = TRUE,
  gg_style = TRUE
)
#> Size at morphometric maturity = 118.6 
#> Confidence intervals = 116 - 121 
#> Rsquare = 0.71

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

p_morph <- plot(
  my_morph_fq,
  xlab = "Carapace width (mm)",
  ylab = "Proportion mature",
  col = c("steelblue", "firebrick"),
  gg_style = TRUE
)
#> Size at morphometric maturity = 118.6 
#> Confidence intervals = 116 - 121 
#> Rsquare = 0.71

names(p_morph)
#> [1] "A"     "B"     "L50"   "ogive"

The individual plots can then be displayed separately.

p_morph$A

p_morph$B

p_morph$L50

p_morph$ogive

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

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.

data(matFish)

head(matFish)
#>   total_length stage_mat
#> 1           12         I
#> 2           12         I
#> 3           13         I
#> 4           14         I
#> 5           14         I
#> 6           14         I

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.

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)
#> formula: Y = 1/[1+exp-(A + B*X)]
#>     Original Bootstrap (Median)
#> A   -8.6047  -8.6497           
#> B   0.356    0.3578            
#> L50 24.1694  24.1841           
#> R2  0.5595   -

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)
#> formula: Y = 1/[1+exp-(A + B*X)]
#>     Bootstrap (Median)
#> A              -8.6026
#> B               0.3565
#> L50            24.1457
#> R2              0.5595

3.2 Base graphics plot for gonadal maturity

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

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")
)

#> Size at gonadal maturity = 24.2 
#> Confidence intervals = 23.8 - 24.6 
#> Rsquare = 0.56

par(oldpar)

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

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

#> Size at gonadal maturity = 24.2 
#> Confidence intervals = 23.8 - 24.6 
#> Rsquare = 0.56

The legend position can be modified with legendPosition.

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

#> Size at gonadal maturity = 24.2 
#> Confidence intervals = 23.8 - 24.6 
#> Rsquare = 0.56

The legend can also be removed.

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

#> Size at gonadal maturity = 24.2 
#> Confidence intervals = 23.8 - 24.6 
#> Rsquare = 0.56

3.3 ggplot2-style plots for gonadal maturity

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

plot(
  my_gonad_fq,
  xlab = "Total length (cm)",
  ylab = "Proportion mature",
  col = c("steelblue", "firebrick"),
  onlyOgive = TRUE,
  gg_style = TRUE
)
#> Size at gonadal maturity = 24.2 
#> Confidence intervals = 23.8 - 24.6 
#> Rsquare = 0.56

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

plot(
  my_gonad_fq,
  xlab = "Total length (cm)",
  ylab = "Proportion mature",
  col = c("steelblue", "firebrick"),
  onlyOgive = TRUE,
  gg_style = TRUE,
  legendPosition = "bottomright"
)
#> Size at gonadal maturity = 24.2 
#> Confidence intervals = 23.8 - 24.6 
#> Rsquare = 0.56

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

p_gonad <- plot(
  my_gonad_fq,
  xlab = "Total length (cm)",
  ylab = "Proportion mature",
  col = c("steelblue", "firebrick"),
  gg_style = TRUE
)
#> Size at gonadal maturity = 24.2 
#> Confidence intervals = 23.8 - 24.6 
#> Rsquare = 0.56

names(p_gonad)
#> [1] "A"     "B"     "L50"   "ogive"

The individual plots can be displayed separately.

p_gonad$A

p_gonad$B

p_gonad$L50

p_gonad$ogive

The returned plots can be customized using ggplot2.

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.

sizeMat: An R Package to Estimate Size at Sexual Maturity

Josymar Torrejon-Magallanes

2026-05-29