---
title: "Length at maturity analysis"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Length at maturity analysis}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 6,
  fig.height = 4
)
library(ggplot2)

# Figures are numbered in the order listed in figs. fig() gives a linked
# reference for the text and fig_cap() a numbered caption. Each captioned
# figure is anchored by its chunk label, which is the link target
figs <- c("mat-female", "mat-proportions", "mat-point", "mat-rug",
          "mat-binwidth", "mat-age")
fig <- function(id) sprintf("[Figure %d](#%s)", match(id, figs), id)
fig_cap <- function(text) {
  sprintf("**Figure %d.** %s", match(knitr::opts_current$get("label"), figs), text)
}
knitr::knit_hooks$set(fig.cap = function(before, options) {
  if (before) sprintf('\n\n<div id="%s"></div>\n\n', options$label)
})
```

```{css, echo = FALSE}
/* Set figure captions apart from the body text. The first selector is used
   by rmarkdown::html_vignette and the second by pkgdown */
.figcaption, p.caption {
  font-size: 0.9em;
  color: #555;
  margin-top: 0.25em;
  margin-bottom: 1.5em;
}
```

The `maturity()` function fits a logistic regression of binary maturity (0 = immature, 1 = mature) as a function of a continuous predictor such as length or age:

$$P(\text{mature}) = \frac{1}{1 + e^{-(a + b \cdot L)}}$$

Key outputs are $L_{50}$ and $L_{95}$ — the lengths (or ages) at which 50% and 95% of individuals are expected to be mature — with bootstrap confidence intervals.

## Basic usage

```{r}
library(mustelus)
data(spottail)

lm50 <- maturity(maturity_stage, length, data = spottail, times = 200)
summary(lm50)
```

## Grouping by sex

```{r}
lm50_sex <- maturity(maturity_stage, length, sex, data = spottail, times = 200)
summary(lm50_sex)
```

## Plotting

`plot()` returns a named list of ggplot objects. The default shows binned observed proportions scaled by sample size, as for females in `r fig("mat-female")`.

```{r mat-female, fig.cap = fig_cap("Proportion mature against total length for female spot-tail sharks, with the fitted logistic curve and its bootstrap 95% confidence interval (grey ribbon). Points are observed proportions in length bins, sized by the number of animals in each bin.")}
p <- plot(lm50_sex)
p$f + xlab("Total length (mm)") + ylab("Proportion mature")
```

### Plot styles

Three styles for raw data display are available, shown for both sexes combined in `r fig("mat-proportions")` to `r fig("mat-rug")`:

```{r mat-proportions, fig.cap = fig_cap(sprintf("Proportion mature against total length for spot-tail sharks of both sexes combined, with points sized by the number of animals in each bin (the default). Curve and ribbon as in %s.", fig("mat-female")))}
# Proportions (default) — point size scaled by n per bin
plot(lm50, raw_data = "proportions")$Unspecified +
  xlab("Total length (mm)") + ylab("Proportion mature")
```

```{r mat-point, fig.cap = fig_cap(sprintf("As %s, with points of uniform size.", fig("mat-proportions")))}
# Point — same bins, uniform point size
plot(lm50, raw_data = "point")$Unspecified +
  xlab("Total length (mm)") + ylab("Proportion mature")
```

```{r mat-rug, fig.cap = fig_cap(sprintf("As %s, with individual animals shown as a rug, mature animals along the top axis and immature animals along the bottom.", fig("mat-proportions")))}
# Rug — mature individuals on top, immature on bottom
plot(lm50, raw_data = "rug")$Unspecified +
  xlab("Total length (mm)") + ylab("Proportion mature")
```

### Adjusting bin width

The `binwidth` argument controls the width of bins in data units (default is range/10). `r fig("mat-binwidth")` uses 50 mm bins:

```{r mat-binwidth, fig.cap = fig_cap(sprintf("As %s, with 50 mm bins.", fig("mat-proportions")))}
plot(lm50, binwidth = 50)$Unspecified +
  xlab("Total length (mm)") + ylab("Proportion mature")
```

## Using age as the predictor

`maturity()` works with any continuous predictor — here using age rather than length:

```{r}
am50 <- maturity(maturity_stage, age_agree, data = spottail, times = 200)
summary(am50)
```

The fitted ogive is shown in `r fig("mat-age")`.

```{r mat-age, fig.cap = fig_cap(sprintf("Proportion mature against age for spot-tail sharks of both sexes combined, in 1 year bins. Curve and ribbon as in %s.", fig("mat-female")))}
plot(am50, binwidth = 1)$Unspecified +
  xlab("Age (years)") + ylab("Proportion mature")
```

## Output structure

Like `len_weight()`, the result is a nested tibble with named list columns per group:

```{r}
# Bootstrap L50 estimates for all individuals
quantile(lm50_sex$mods[["f"]]$fitted.values)

# Predicted maturity curve for females
head(lm50_sex$preds[["f"]])
```