---
title: "Gene-Drop Iteration Convergence"
subtitle: "nprcgenekeepr: an R Package for the Genetic Management of Colonies"
author: "R. Mark Sharp, Ph.D."
date: "2026-06-25"
output:
  html_document:
    df_print: paged
vignette: >
  %\VignetteEngine{knitr::rmarkdown_notangle}
  %\VignetteIndexEntry{Gene-Drop Iteration Convergence}
  %\usepackage[UTF-8]{inputenc}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
library(nprcgenekeepr)
library(knitr)
library(kableExtra)
```

## How many gene-drop iterations does a pedigree need?

The Genetic Value Analysis ranks animals by mean kinship and **genome
uniqueness**. Mean kinship is deterministic, but genome uniqueness (`gu`) is
*estimated* by a gene-drop simulation: alleles are dropped down the pedigree a
chosen number of times, and `gu` is the average result. Like any estimate from a
finite number of draws, it carries sampling noise, so the precision of `gu` --
and, more importantly, the **selection order** it produces -- depends on how many
iterations you run.

There is no single iteration count that is right for every colony. A pedigree
with many animals whose true genome uniqueness sits near a ranking boundary needs
more iterations before the chosen set of animals stops changing run to run; a
pedigree whose ranked animals are well separated settles almost immediately.
`gvaConvergence()` answers the question *for the pedigree in hand*: it reports how
reproducible the selection order is at a range of iteration counts, and the
smallest count at which it has settled.

### What it measures

Because the gene-drop iteration columns are independent replicates, the whole
convergence picture is recoverable from a **single** gene drop. `gvaConvergence()`
runs one gene drop at a budget `nMax`, then for each candidate iteration count `N`
it splits the columns into two disjoint halves of `N` columns each -- two genuinely
independent `N`-iteration runs -- ranks each half through the same ordering pipeline
the report uses, and compares the two orderings on two criteria:

- **top-`k` overlap** -- do the two runs choose the same top animals? (default
  `k = 20`, threshold `oMin = 0.90`)
- **Kendall rank agreement** -- do the chosen animals come out in the same order?
  (threshold `rhoMin = 0.95`)

The run is **reproducible at `N`** when both hold, and the recommended iteration
count is the smallest `N` meeting both. The "Undetermined" animals (both
parents unknown, no recorded origin) are a policy constant with rank `NA`; they are
excluded from the order and reported separately as `nUndetermined`.

This is distinct from seed reproducibility. A fixed seed already makes `gu`
bit-identical run to run -- that is reproducibility of the *process*.
`gvaConvergence()` reports the sampling reproducibility of the *estimate*: would a
fresh, independent run lead to the same animals being chosen, in the same order?

## A pedigree where the iteration count matters

No bundled pedigree exercises the diagnostic well -- after the `gu = 0`
de-inflation the shipped pedigrees have no `gu` signal left to rank on, so their selection
order is settled at the smallest iteration count. To show a pedigree where the
count *does* matter, we build a small half-sib web in which the founders (the
sources of the private alleles) are excluded from the analyzed population, so their
alleles survive among the offspring only through descendants. Overlapping
sire/dam mating windows give the offspring distinct true genome-uniqueness values
that straddle the ranking boundary, so at low iteration counts the gene-drop
estimate randomly crosses the boundary and the selection order churns.

```{r fixture-build}
## A deterministic dense-mid-range pedigree: 14 founder sires, each mated to a
## wrapping window of 5 of 15 founder dams (one offspring per pair). Windows
## overlap, so dam fan sizes vary and the 70 offspring get distinct genome
## uniqueness straddling the 10% ranking cutoff. Founders are excluded from the
## analyzed population (`pop`); the offspring are the probands.
makeConvergenceFixture <- function() {
  w <- rep(5L, 14L)          # sire window widths
  b <- 15L                   # number of founder dams
  a <- length(w)
  sids <- sprintf("S%03d", seq_len(a))
  dids <- sprintf("D%03d", seq_len(b))
  id <- c(sids, dids)
  sire <- rep(NA_character_, a + b)
  dam <- rep(NA_character_, a + b)
  sex <- c(rep("M", a), rep("F", b))
  off <- character(0L)
  ocount <- 0L
  start <- 1L
  for (i in seq_len(a)) {
    for (j in seq_len(w[i])) {
      dj <- ((start + j - 2L) %% b) + 1L
      ocount <- ocount + 1L
      o <- sprintf("O%04d", ocount)
      id <- c(id, o)
      sire <- c(sire, sids[i])
      dam <- c(dam, dids[dj])
      sex <- c(sex, if (ocount %% 2L == 0L) "F" else "M")
      off <- c(off, o)
    }
    start <- start + max(1L, w[i] - 1L)
  }
  ped <- data.frame(id = id, sire = sire, dam = dam, sex = sex,
                    stringsAsFactors = FALSE)
  ped$gen <- findGeneration(ped$id, ped$sire, ped$dam)
  list(ped = ped, pop = off)
}

fx <- makeConvergenceFixture()
nrow(fx$ped)      # total animals (founders + offspring)
length(fx$pop)    # offspring (the analyzed population)
```

We assess iteration counts from 25 up to 800 (a fixed `seed` makes the curve
reproducible):

```{r fixture-convergence}
conv <- gvaConvergence(
  fx$ped, pop = fx$pop, nMax = 1600L,
  grid = c(25L, 50L, 100L, 200L, 400L, 800L), seed = 11L
)

kable(
  conv$convergence,
  digits = c(0L, 3L, 3L),
  col.names = c("Iterations (N)", "Top-20 overlap", "Kendall agreement"),
  caption = "Selection-order reproducibility vs. iteration count (half-sib web)."
) |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

conv$recommendedIter   # smallest N meeting both criteria
conv$converged
conv$nRankable
conv$nUndetermined
```

At the smallest counts the two independent half-runs disagree on both which
animals are chosen and their order; both measures climb as `N` grows, and the
selection order has settled by the higher counts and stays settled there.
`recommendedIter` (printed below the table) reports the smallest count meeting
both criteria. For this pedigree the default of 1000 iterations is in the right
range, and `gvaConvergence()` is how you would confirm that rather than guess it.

## A pedigree that converges immediately

Contrast the bundled `qcPed`. After the `gu = 0` de-inflation none of its
ranked animals carry a non-zero genome uniqueness, so the selection order is
driven by deterministic mean kinship and does not move with the iteration count.

```{r qcped-convergence}
convQc <- gvaConvergence(
  nprcgenekeepr::qcPed, nMax = 400L,
  grid = c(25L, 50L, 100L, 200L), seed = 11L
)

kable(
  convQc$convergence,
  digits = c(0L, 3L, 3L),
  col.names = c("Iterations (N)", "Top-20 overlap", "Kendall agreement"),
  caption = "Selection-order reproducibility vs. iteration count (qcPed)."
) |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

convQc$recommendedIter   # converges at the grid floor
convQc$nRankable
convQc$nUndetermined     # the excluded Undetermined set
```

Overlap and agreement are at their maximum from the smallest count, so the
recommended count is the floor of the grid: this pedigree needs almost no
iterations for the *order* to be reproducible. The `nUndetermined` count reports
how many animals were set aside without a rank.

## Choosing an iteration count

Run `gvaConvergence()` on your own pedigree and read `recommendedIter`. If
`converged` is `FALSE`, no count in the grid settled the order -- raise `nMax` and
extend the `grid` (each candidate `N` needs `2 * N <= nMax`), or accept that this
pedigree is intrinsically near a ranking boundary and interpret the ranking with
care. The thresholds `k`, `oMin`, and `rhoMin` are arguments, so you can make the
definition of "reproducible" stricter or looser for your purpose.

Keep two ideas separate. The per-animal sampling standard error (`guSE`, reported
beside `gu` in `reportGV()`) tells you how *precise the genome uniqueness number*
is; it shrinks like one over the square root of the iteration count. The
*selection order* -- what `gvaConvergence()` measures -- is what actually decides
which animals are chosen for breeding. A small `guSE` does not by itself mean the
order has settled, which is exactly why the convergence check exists.

See also `reportGV()` (the analysis whose iterations this advises on),
`calcGUSE()` (the per-animal standard error), and the *Colony Manager Tutorial*
vignette.