The coalescent and gene trees

population genetics

ecology tutorial

ggplot2

Simulate the coalescent in R, check the tree height and total branch length against theory, and see why a single gene tree says almost nothing about history.

Author

Tidy Ecology

Published

2026-07-11

Modified

2026-10-14

Corrected 14 October 2026: the example gene tree placed its tips in sampling order, so several branches crossed junction bars and two of the three long lineages were drawn almost on top of each other; the tips now follow the tree’s branching order, and the tree itself is unchanged. The text called the average per-tree share of the two-lineage interval slightly lower than the ratio of expectations, although it falls below one half while the ratio is above it; the word slightly has gone.

Corrected 7 October 2026: the post said the tree height can never exceed two coalescent units however many individuals are sequenced; that bound holds for the expected height, while a single tree can be taller, and the sentence now says so. The expected standard deviation of the pairwise coalescence time in the Wright-Fisher check was given as the mean of the geometric distribution and is now that distribution’s own standard deviation, a little below the mean.

Corrected 27 September 2026: the caption of the example tree said its final two-lineage interval takes up much of the height; in the tree shown that interval is short and three lineages fill most of the depth, and the caption now says so. The text before it now calls the half-the-depth share an average over trees rather than the shape of every tree, and the preview image, whose subtitle said half the depth is the last interval, now says three lineages fill most of this tree.

Everything in this cluster so far has run forwards: start with a population, let drift and mutation act, see what the sample looks like. The coalescent runs the same process backwards. Take the sample, trace ancestral lineages until they meet, and you get a gene tree whose shape carries the demographic history that produced it.

The backward view is not just an alternative; it is what makes computation cheap. A Wright-Fisher simulation has to track a whole population to say something about ten sampled gene copies. The coalescent tracks the ten lineages and nothing else, and one line of rexp() does it.

Waiting for a coalescence

With k ancestral lineages in a population of 2N gene copies, any particular pair shares a parent in the previous generation with probability 1/(2N), and there are k(k-1)/2 pairs. Measure time in units of 2N generations and the waiting time until the next coalescence is exponential with rate k(k-1)/2. That is the whole model.

library(ggplot2)

te_pal <- list(forest = "#275139", green = "#2f8f63", sage = "#93a87f",
               clay = "#b5534e", gold = "#cda23f", line = "#dad9ca",
               ink = "#16241d")

theme_te <- function() {
  theme_minimal(base_size = 12) +
    theme(panel.grid.minor = element_blank(),
          panel.grid.major = element_line(colour = "#e7e6dc"),
          plot.background = element_rect(fill = "#f5f4ee", colour = NA),
          panel.background = element_rect(fill = "#f5f4ee", colour = NA),
          plot.title = element_text(face = "bold", colour = te_pal$ink),
          axis.title = element_text(colour = "#2c3a31"))
}

coal_times <- function(n) {
  k <- n:2
  rexp(length(k), rate = k * (k - 1) / 2)
}

tmrca <- function(n) sum(coal_times(n))

total_len <- function(n) {
  k <- n:2
  sum(rexp(length(k), rate = k * (k - 1) / 2) * k)
}

set.seed(477)
reps <- 20000
n_grid <- c(2, 5, 10, 20, 50)
tab <- t(sapply(n_grid, function(n) {
  tt <- replicate(reps, tmrca(n))
  tl <- replicate(reps, total_len(n))
  c(n = n, tmrca_sim = mean(tt), tmrca_theory = 2 * (1 - 1 / n),
    tmrca_sd = sd(tt), total_sim = mean(tl),
    total_theory = 2 * sum(1 / seq_len(n - 1)))
}))
print(round(tab, 3))

      n tmrca_sim tmrca_theory tmrca_sd total_sim total_theory
[1,]  2     1.006         1.00    0.990     2.035        2.000
[2,]  5     1.606         1.60    1.076     4.176        4.167
[3,] 10     1.806         1.80    1.092     5.624        5.658
[4,] 20     1.891         1.90    1.077     7.084        7.095
[5,] 50     1.961         1.96    1.073     8.948        8.958

Two expectations fall out of the exponential waiting times, and the simulation reproduces both. The time to the most recent common ancestor is 2(1 - 1/n) in units of 2N generations, so 1.006 against 1.000 for a pair and 1.961 against 1.960 for fifty lineages. The total branch length is 2 * sum(1/k) for k from 1 to n-1, giving 8.948 against 8.958 at n = 50.

Look at how differently those two quantities grow. Going from two samples to fifty roughly doubles the expected tree height, from 1.0 to 1.96, and that expectation can never exceed 2 no matter how many individuals are sequenced, although a single tree can be taller. The total branch length grows without bound, but only logarithmically: 2.0, 4.167, 5.658, 7.095, 8.958. Mutations accumulate in proportion to total branch length, so this is the coalescent statement of the sampling result from the previous tutorial: the tenth individual adds a lot, the fiftieth adds very little, and neither adds much depth.

One tree is one draw

The standard deviation column is the part that gets forgotten. At n = 10 the mean TMRCA is 1.806 and the standard deviation is 1.092.

set.seed(478)
n_show <- 10
parts <- t(replicate(5000, {
  ts <- coal_times(n_show)
  c(t2_share = ts[length(ts)] / sum(ts), tmrca = sum(ts))
}))

cat("mean share of the tree height spent with two lineages:",
    round(mean(parts[, "t2_share"]), 4),
    " ratio of expectations 1/(2(1 - 1/n)):",
    round(1 / (2 * (1 - 1 / n_show)), 4), "\n")

mean share of the tree height spent with two lineages: 0.4683  ratio of expectations 1/(2(1 - 1/n)): 0.5556

cat("TMRCA coefficient of variation at n = 10:",
    round(sd(parts[, "tmrca"]) / mean(parts[, "tmrca"]), 4), "\n")

TMRCA coefficient of variation at n = 10: 0.604

A coefficient of variation of 0.604 means the tree height for one locus is essentially a single draw from a very broad distribution. Two loci with the same demographic history can easily differ threefold in depth, and no amount of extra sampling within a locus fixes that, because all the sampled copies share the same tree. This is the deepest structural reason why population genetic inference needs many independent loci rather than deep sampling of a few.

The other number describes the typical shape of a gene tree. The interval during which only two lineages remain has expectation 1, out of a total expectation of 1.8, so 0.5556 of the tree height sits in that last interval. Averaging the per-tree share rather than the expectations gives 0.4683, lower because the ratio is not linear. Either way, on average roughly half the depth of a gene tree is one long stretch with two lineages, and the coalescences that separate the ten samples are crowded into the recent end. Any single tree can sit a long way from that average, as the one below shows.

set.seed(480)

draw_tree <- function(n) {
  # simulate the merges first: one exponential wait and one random pair per step
  active <- as.list(seq_len(n))   # the tips below each remaining lineage
  height <- 0
  merges <- list()
  while (length(active) > 1) {
    k <- length(active)
    height <- height + rexp(1, rate = k * (k - 1) / 2)
    pick <- sample(k, 2)
    merges[[length(merges) + 1]] <- list(pick = pick, h = height)
    keep <- setdiff(seq_len(k), pick)
    active <- c(active[keep], list(unlist(active[pick])))
  }
  # each clade's tips sit together in the root's tip list, so placing the
  # tips in that order keeps every branch clear of the junction bars
  pos <- match(seq_len(n), active[[1]])
  base <- rep(0, n)
  segs <- data.frame()
  for (m in merges) {
    pick <- m$pick
    segs <- rbind(segs,
                  data.frame(x = pos[pick], xend = pos[pick], y = base[pick], yend = m$h),
                  data.frame(x = pos[pick[1]], xend = pos[pick[2]], y = m$h, yend = m$h))
    keep <- setdiff(seq_along(pos), pick)
    pos <- c(pos[keep], mean(pos[pick]))
    base <- c(base[keep], m$h)
  }
  segs
}
segs <- draw_tree(10)
cat("this tree: height", round(max(segs$yend), 3), "\n")

this tree: height 2.765

ggplot(segs) +
  geom_segment(aes(x = x, xend = xend, y = y, yend = yend),
               colour = te_pal$forest, linewidth = 0.7) +
  labs(title = "A gene tree for ten sampled copies",
       subtitle = "coalescent units of 2N generations; time runs upward into the past",
       x = "sampled lineages", y = "time before the present") +
  theme_te() +
  theme(axis.text.x = element_blank(), panel.grid.major.x = element_blank())

A gene tree growing upward from ten tips; seven junctions are crowded near the tips, three long branches run up most of the height, and only a short final stretch with two lineages remains below the root. — Figure 1: One realised coalescent tree for ten sampled lineages, drawn on the coalescent timescale. In this draw three lineages persist for most of the height and the final two-lineage interval is short: the half-the-depth figure is an average over trees, and one tree is one draw.

theory_df <- data.frame(n = 2:60)
theory_df$tmrca <- 2 * (1 - 1 / theory_df$n)
theory_df$total <- sapply(theory_df$n, function(n) 2 * sum(1 / seq_len(n - 1)))
long_theory <- rbind(
  data.frame(n = theory_df$n, value = theory_df$tmrca, what = "tree height"),
  data.frame(n = theory_df$n, value = theory_df$total,
             what = "total branch length"))
sim_df <- as.data.frame(tab)
long_sim <- rbind(
  data.frame(n = sim_df$n, value = sim_df$tmrca_sim, what = "tree height"),
  data.frame(n = sim_df$n, value = sim_df$total_sim,
             what = "total branch length"))

ggplot(long_theory, aes(n, value, colour = what)) +
  geom_line(linewidth = 0.9) +
  geom_point(data = long_sim, size = 2.4) +
  scale_x_log10(breaks = c(2, 5, 10, 20, 50)) +
  scale_colour_manual(values = c(`tree height` = te_pal$clay,
                                 `total branch length` = te_pal$forest),
                      name = NULL) +
  labs(title = "Depth saturates, branch length does not",
       subtitle = "lines: coalescent expectations, points: 20000 simulated trees per sample size",
       x = "sampled lineages (log scale)",
       y = "coalescent units of 2N generations") +
  theme_te()

Two rising curves on a log-x axis; tree height flattens below two while total branch length keeps climbing logarithmically, with simulated points on both theory lines. — Figure 2: Tree height and total branch length against sample size, simulated against the coalescent expectations.

Does it match a real population?

The coalescent is an approximation to the Wright-Fisher model, valid when the sample is small relative to the population. It is worth checking that against the forward simulation the earlier tutorials used, rather than taking it on faith. Two gene copies in a population of N diploids should coalesce after a geometric number of generations with mean 2N.

set.seed(479)
wf_pair <- function(n_dip, max_gen = 20000) {
  for (g in seq_len(max_gen)) {
    pa <- sample(2 * n_dip, 1)
    pb <- sample(2 * n_dip, 1)
    if (pa == pb) return(g)
  }
  NA
}

n_dip <- 50
pair_times <- replicate(3000, wf_pair(n_dip))
cat("Wright-Fisher with N =", n_dip,
    ": mean pairwise coalescence time", round(mean(pair_times), 1),
    " expected 2N", 2 * n_dip, "\n")

Wright-Fisher with N = 50 : mean pairwise coalescence time 98.2  expected 2N 100

cat("standard deviation", round(sd(pair_times), 1),
    " expected (geometric)", round(sqrt(1 - 1 / (2 * n_dip)) * 2 * n_dip, 1),
    " in coalescent units", round(mean(pair_times) / (2 * n_dip), 3), "\n")

standard deviation 96.6  expected (geometric) 99.5  in coalescent units 0.982

Mean 98.2 generations against an expected 100, with a standard deviation of 96.6 against an expected 99.5: the geometric distribution has a standard deviation almost equal to its mean, which is where the exponential approximation comes from and why coalescence times are so variable. Dividing by 2N puts the forward simulation on the coalescent timescale at 0.982, and the two models agree.

That correspondence is what makes Ne usable in the coalescent. The backward process does not know about sex ratios or reproductive variance; it only knows a rate. Everything the previous tutorials did to Ne shows up here as a change of timescale, which is why coalescent methods estimate an effective size rather than a census.

Where to go next

A tree by itself is not data. What gets observed is mutations, and mutations fall on branches in proportion to their length, which turns the tree into a pattern of variable sites. The next tutorial adds mutation to these trees and derives the two classical diversity estimators, one counting differences between sequences and one counting segregating sites, then measures which is more precise and why.

References

Kingman JFC 1982 Stochastic Processes and their Applications 13(3):235-248 (10.1016/0304-4149(82)90011-4)

Hudson RR 1990 Oxford Surveys in Evolutionary Biology 7:1-44 (no DOI)

Rosenberg NA, Nordborg M 2002 Nature Reviews Genetics 3(5):380-390 (10.1038/nrg795)

Wakeley J 2009 Coalescent Theory: An Introduction, Macmillan, ISBN 978-0974707754

--- title: "The coalescent and gene trees" description: "Simulate the coalescent in R, check the tree height and total branch length against theory, and see why a single gene tree says almost nothing about history." date: "2026-07-11 11:00" categories: [R, population genetics, ecology tutorial, ggplot2] image: thumbnail.png image-alt: "A simulated gene tree for ten sampled lineages, with seven coalescences packed close to the present, three lineages left for most of the depth, and only a short final interval where two lineages remain." date-modified: "2026-10-14" --- *Corrected 14 October 2026: the example gene tree placed its tips in sampling order, so several branches crossed junction bars and two of the three long lineages were drawn almost on top of each other; the tips now follow the tree's branching order, and the tree itself is unchanged. The text called the average per-tree share of the two-lineage interval slightly lower than the ratio of expectations, although it falls below one half while the ratio is above it; the word slightly has gone.* *Corrected 7 October 2026: the post said the tree height can never exceed two coalescent units however many individuals are sequenced; that bound holds for the expected height, while a single tree can be taller, and the sentence now says so. The expected standard deviation of the pairwise coalescence time in the Wright-Fisher check was given as the mean of the geometric distribution and is now that distribution's own standard deviation, a little below the mean.* *Corrected 27 September 2026: the caption of the example tree said its final two-lineage interval takes up much of the height; in the tree shown that interval is short and three lineages fill most of the depth, and the caption now says so. The text before it now calls the half-the-depth share an average over trees rather than the shape of every tree, and the preview image, whose subtitle said half the depth is the last interval, now says three lineages fill most of this tree.* Everything in this cluster so far has run forwards: start with a population, let drift and mutation act, see what the sample looks like. The coalescent runs the same process backwards. Take the sample, trace ancestral lineages until they meet, and you get a gene tree whose shape carries the demographic history that produced it. The backward view is not just an alternative; it is what makes computation cheap. A Wright-Fisher simulation has to track a whole population to say something about ten sampled gene copies. The coalescent tracks the ten lineages and nothing else, and one line of `rexp()` does it. ## Waiting for a coalescence With `k` ancestral lineages in a population of `2N` gene copies, any particular pair shares a parent in the previous generation with probability `1/(2N)`, and there are `k(k-1)/2` pairs. Measure time in units of `2N` generations and the waiting time until the next coalescence is exponential with rate `k(k-1)/2`. That is the whole model. ```{r setup} #| message: false #| warning: false library(ggplot2) te_pal <- list(forest = "#275139", green = "#2f8f63", sage = "#93a87f", clay = "#b5534e", gold = "#cda23f", line = "#dad9ca", ink = "#16241d") theme_te <- function() { theme_minimal(base_size = 12) + theme(panel.grid.minor = element_blank(), panel.grid.major = element_line(colour = "#e7e6dc"), plot.background = element_rect(fill = "#f5f4ee", colour = NA), panel.background = element_rect(fill = "#f5f4ee", colour = NA), plot.title = element_text(face = "bold", colour = te_pal$ink), axis.title = element_text(colour = "#2c3a31")) } ``` ```{r coalescent} coal_times <- function(n) { k <- n:2 rexp(length(k), rate = k * (k - 1) / 2) } tmrca <- function(n) sum(coal_times(n)) total_len <- function(n) { k <- n:2 sum(rexp(length(k), rate = k * (k - 1) / 2) * k) } set.seed(477) reps <- 20000 n_grid <- c(2, 5, 10, 20, 50) tab <- t(sapply(n_grid, function(n) { tt <- replicate(reps, tmrca(n)) tl <- replicate(reps, total_len(n)) c(n = n, tmrca_sim = mean(tt), tmrca_theory = 2 * (1 - 1 / n), tmrca_sd = sd(tt), total_sim = mean(tl), total_theory = 2 * sum(1 / seq_len(n - 1))) })) print(round(tab, 3)) ``` Two expectations fall out of the exponential waiting times, and the simulation reproduces both. The time to the most recent common ancestor is `2(1 - 1/n)` in units of `2N` generations, so `r sprintf("%.3f", tab[1, "tmrca_sim"])` against `r sprintf("%.3f", tab[1, "tmrca_theory"])` for a pair and `r sprintf("%.3f", tab[5, "tmrca_sim"])` against `r sprintf("%.3f", tab[5, "tmrca_theory"])` for fifty lineages. The total branch length is `2 * sum(1/k)` for `k` from 1 to `n-1`, giving `r sprintf("%.3f", tab[5, "total_sim"])` against `r sprintf("%.3f", tab[5, "total_theory"])` at `n` = 50. Look at how differently those two quantities grow. Going from two samples to fifty roughly doubles the expected tree height, from `r sprintf("%.1f", tab[1, "tmrca_theory"])` to `r sprintf("%.2f", tab[5, "tmrca_theory"])`, and that expectation can never exceed 2 no matter how many individuals are sequenced, although a single tree can be taller. The total branch length grows without bound, but only logarithmically: `r sprintf("%.1f", tab[1, "total_theory"])`, `r sprintf("%.3f", tab[2, "total_theory"])`, `r sprintf("%.3f", tab[3, "total_theory"])`, `r sprintf("%.3f", tab[4, "total_theory"])`, `r sprintf("%.3f", tab[5, "total_theory"])`. Mutations accumulate in proportion to total branch length, so this is the coalescent statement of the sampling result from the previous tutorial: the tenth individual adds a lot, the fiftieth adds very little, and neither adds much depth. ## One tree is one draw The standard deviation column is the part that gets forgotten. At `n` = 10 the mean `TMRCA` is `r sprintf("%.3f", tab[3, "tmrca_sim"])` and the standard deviation is `r sprintf("%.3f", tab[3, "tmrca_sd"])`. ```{r variance} set.seed(478) n_show <- 10 parts <- t(replicate(5000, { ts <- coal_times(n_show) c(t2_share = ts[length(ts)] / sum(ts), tmrca = sum(ts)) })) cat("mean share of the tree height spent with two lineages:", round(mean(parts[, "t2_share"]), 4), " ratio of expectations 1/(2(1 - 1/n)):", round(1 / (2 * (1 - 1 / n_show)), 4), "\n") cat("TMRCA coefficient of variation at n = 10:", round(sd(parts[, "tmrca"]) / mean(parts[, "tmrca"]), 4), "\n") ``` A coefficient of variation of `r sprintf("%.3f", sd(parts[, "tmrca"]) / mean(parts[, "tmrca"]))` means the tree height for one locus is essentially a single draw from a very broad distribution. Two loci with the same demographic history can easily differ threefold in depth, and no amount of extra sampling within a locus fixes that, because all the sampled copies share the same tree. This is the deepest structural reason why population genetic inference needs many independent loci rather than deep sampling of a few. The other number describes the typical shape of a gene tree. The interval during which only two lineages remain has expectation 1, out of a total expectation of `r sprintf("%.1f", 2 * (1 - 1 / n_show))`, so `r sprintf("%.4f", 1 / (2 * (1 - 1 / n_show)))` of the tree height sits in that last interval. Averaging the per-tree share rather than the expectations gives `r sprintf("%.4f", mean(parts[, "t2_share"]))`, lower because the ratio is not linear. Either way, on average roughly half the depth of a gene tree is one long stretch with two lineages, and the coalescences that separate the ten samples are crowded into the recent end. Any single tree can sit a long way from that average, as the one below shows. ```{r fig-tree} #| fig-width: 7 #| fig-height: 4.4 #| fig-cap: "One realised coalescent tree for ten sampled lineages, drawn on the coalescent timescale. In this draw three lineages persist for most of the height and the final two-lineage interval is short: the half-the-depth figure is an average over trees, and one tree is one draw." #| fig-alt: "A gene tree growing upward from ten tips; seven junctions are crowded near the tips, three long branches run up most of the height, and only a short final stretch with two lineages remains below the root." set.seed(480) draw_tree <- function(n) { # simulate the merges first: one exponential wait and one random pair per step active <- as.list(seq_len(n)) # the tips below each remaining lineage height <- 0 merges <- list() while (length(active) > 1) { k <- length(active) height <- height + rexp(1, rate = k * (k - 1) / 2) pick <- sample(k, 2) merges[[length(merges) + 1]] <- list(pick = pick, h = height) keep <- setdiff(seq_len(k), pick) active <- c(active[keep], list(unlist(active[pick]))) } # each clade's tips sit together in the root's tip list, so placing the # tips in that order keeps every branch clear of the junction bars pos <- match(seq_len(n), active[[1]]) base <- rep(0, n) segs <- data.frame() for (m in merges) { pick <- m$pick segs <- rbind(segs, data.frame(x = pos[pick], xend = pos[pick], y = base[pick], yend = m$h), data.frame(x = pos[pick[1]], xend = pos[pick[2]], y = m$h, yend = m$h)) keep <- setdiff(seq_along(pos), pick) pos <- c(pos[keep], mean(pos[pick])) base <- c(base[keep], m$h) } segs } segs <- draw_tree(10) cat("this tree: height", round(max(segs$yend), 3), "\n") ggplot(segs) + geom_segment(aes(x = x, xend = xend, y = y, yend = yend), colour = te_pal$forest, linewidth = 0.7) + labs(title = "A gene tree for ten sampled copies", subtitle = "coalescent units of 2N generations; time runs upward into the past", x = "sampled lineages", y = "time before the present") + theme_te() + theme(axis.text.x = element_blank(), panel.grid.major.x = element_blank()) ``` ```{r fig-scaling} #| fig-width: 7 #| fig-height: 4 #| fig-cap: "Tree height and total branch length against sample size, simulated against the coalescent expectations." #| fig-alt: "Two rising curves on a log-x axis; tree height flattens below two while total branch length keeps climbing logarithmically, with simulated points on both theory lines." theory_df <- data.frame(n = 2:60) theory_df$tmrca <- 2 * (1 - 1 / theory_df$n) theory_df$total <- sapply(theory_df$n, function(n) 2 * sum(1 / seq_len(n - 1))) long_theory <- rbind( data.frame(n = theory_df$n, value = theory_df$tmrca, what = "tree height"), data.frame(n = theory_df$n, value = theory_df$total, what = "total branch length")) sim_df <- as.data.frame(tab) long_sim <- rbind( data.frame(n = sim_df$n, value = sim_df$tmrca_sim, what = "tree height"), data.frame(n = sim_df$n, value = sim_df$total_sim, what = "total branch length")) ggplot(long_theory, aes(n, value, colour = what)) + geom_line(linewidth = 0.9) + geom_point(data = long_sim, size = 2.4) + scale_x_log10(breaks = c(2, 5, 10, 20, 50)) + scale_colour_manual(values = c(`tree height` = te_pal$clay, `total branch length` = te_pal$forest), name = NULL) + labs(title = "Depth saturates, branch length does not", subtitle = "lines: coalescent expectations, points: 20000 simulated trees per sample size", x = "sampled lineages (log scale)", y = "coalescent units of 2N generations") + theme_te() ``` ## Does it match a real population? The coalescent is an approximation to the Wright-Fisher model, valid when the sample is small relative to the population. It is worth checking that against the forward simulation the earlier tutorials used, rather than taking it on faith. Two gene copies in a population of `N` diploids should coalesce after a geometric number of generations with mean `2N`. ```{r wf-check} set.seed(479) wf_pair <- function(n_dip, max_gen = 20000) { for (g in seq_len(max_gen)) { pa <- sample(2 * n_dip, 1) pb <- sample(2 * n_dip, 1) if (pa == pb) return(g) } NA } n_dip <- 50 pair_times <- replicate(3000, wf_pair(n_dip)) cat("Wright-Fisher with N =", n_dip, ": mean pairwise coalescence time", round(mean(pair_times), 1), " expected 2N", 2 * n_dip, "\n") cat("standard deviation", round(sd(pair_times), 1), " expected (geometric)", round(sqrt(1 - 1 / (2 * n_dip)) * 2 * n_dip, 1), " in coalescent units", round(mean(pair_times) / (2 * n_dip), 3), "\n") ``` Mean `r sprintf("%.1f", mean(pair_times))` generations against an expected 100, with a standard deviation of `r sprintf("%.1f", sd(pair_times))` against an expected `r sprintf("%.1f", sqrt(1 - 1 / (2 * n_dip)) * 2 * n_dip)`: the geometric distribution has a standard deviation almost equal to its mean, which is where the exponential approximation comes from and why coalescence times are so variable. Dividing by `2N` puts the forward simulation on the coalescent timescale at `r sprintf("%.3f", mean(pair_times) / (2 * n_dip))`, and the two models agree. That correspondence is what makes `Ne` usable in the coalescent. The backward process does not know about sex ratios or reproductive variance; it only knows a rate. Everything the previous tutorials did to `Ne` shows up here as a change of timescale, which is why coalescent methods estimate an effective size rather than a census. ## Where to go next A tree by itself is not data. What gets observed is mutations, and mutations fall on branches in proportion to their length, which turns the tree into a pattern of variable sites. The next tutorial adds mutation to these trees and derives the two classical diversity estimators, one counting differences between sequences and one counting segregating sites, then measures which is more precise and why. ## References Kingman JFC 1982 Stochastic Processes and their Applications 13(3):235-248 (10.1016/0304-4149(82)90011-4) Hudson RR 1990 Oxford Surveys in Evolutionary Biology 7:1-44 (no DOI) Rosenberg NA, Nordborg M 2002 Nature Reviews Genetics 3(5):380-390 (10.1038/nrg795) Wakeley J 2009 Coalescent Theory: An Introduction, Macmillan, ISBN 978-0974707754 ## Related tutorials - [Nucleotide diversity and Watterson theta](../nucleotide-diversity-and-watterson-theta/) - [Tajima's D and the frequency spectrum](../tajimas-d-and-the-frequency-spectrum/) - [Checking a coalescent analysis](../checking-a-coalescent-analysis/) - [Effective population size](../effective-population-size/)