---
title: "Fecundity analysis"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Fecundity 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("fec-length", "fec-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;
}
```

## Background

Fecundity in chondrichthyans — the number of *in utero* embryos carried by a
gravid female — typically increases with maternal body size. `fecundity()` fits
the simplest description of that relationship, an ordinary linear regression:

$$F_i = a + b \cdot L_i + \epsilon, \quad \epsilon \sim N(0, \sigma^2)$$

where $F_i$ is the number of embryos carried by female $i$ and $L_i$ is her
length (or age). The slope $b$ is the increase in litter size per unit length,
and is usually the quantity of interest.

## Non-gravid animals are excluded

A fecundity count of zero means an animal was not gravid, not that its fecundity
was zero. Including such records would bias the intercept downwards and inflate
the slope. `fecundity()` therefore drops records with zero or negative counts
before fitting and reports how many were removed. Predictions span the length or
age range of gravid animals only.

## Basic usage

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

fec <- fecundity(emb, length, data = spottail)
summary(fec)
```

The summary reports the intercept $a$ and slope $b$ with their standard errors,
the coefficient of determination, the p-value for the slope, the residual
standard deviation $\sigma$, sample size, mean and range of observed fecundity,
and the predictor range over which the model applies.

## Plotting

`plot()` returns a ggplot object, which can be customised with standard ggplot2
calls. In `r fig("fec-length")` the solid line is the fitted mean, the dashed
ribbon the 95% confidence interval, and the dotted ribbon the 95% prediction
interval — the presentation used by Walker (2005).

```{r fec-length, fig.cap = fig_cap("Number of embryos against maternal total length for gravid spot-tail sharks. The solid line is the fitted mean, the dashed ribbon the 95% confidence interval and the dotted ribbon the 95% prediction interval.")}
plot(fec) + xlab("Total length (mm)") + ylab("Number of embryos")
```

As in `len_weight()`, the two interval types answer different questions. The
confidence interval describes uncertainty in *mean* litter size at a given
length and narrows as data accumulate; the prediction interval describes the
spread of litter sizes you would expect to observe in *individual* females, and
does not narrow appreciably. For predicting the reproductive output of a female
of known size, the prediction interval is the relevant one.

## Using age as the predictor

Any continuous predictor works — here maternal age rather than length:

```{r}
fec_age <- fecundity(emb, age_agree, data = spottail)
summary(fec_age)
```

The fitted relationship is shown in `r fig("fec-age")`.

```{r fec-age, fig.cap = fig_cap(sprintf("Number of embryos against maternal age for gravid spot-tail sharks, with lines as in %s.", fig("fec-length")))}
plot(fec_age) + xlab("Age (years)") + ylab("Number of embryos")
```

## Output structure

Unlike `len_weight()` and `maturity()`, `fecundity()` fits a single model and
takes no grouping variable. The result is a one row tibble whose list columns
hold the components of that fit:

```{r}
# Fitted lm object
fec$mods[[1]]
```

```{r}
# Predicted values
head(fec$preds[[1]])
```

The `preds` data frame contains six columns: `x`, `fec` (predicted mean),
`clower`/`cupper` (95% confidence interval), and `plower`/`pupper` (95%
prediction interval).

## Reference

Walker, T.I. (2005) Reproduction in fisheries science. In: Hamlett, W.C. (ed)
*Reproductive Biology and Phylogeny of Chondrichthyes*.