Multilevel selection and group structure

evolutionary ecology
evolutionary game theory
R
ecology tutorial
Split the Price equation into between-group and within-group terms in R, and build a case where cooperators lose inside every group yet rise in the population.
Author

Tidy Ecology

Published

2026-08-09

The Price equation splits change in a trait into a selection term and a transmission term. It can be split a second time. If the population is divided into groups, the covariance between fitness and trait value can be written as a covariance between group means plus an average of the covariances inside groups. That is the law of total covariance, applied to evolution, and it gives an exact accounting identity:

\[\bar{w}\,\Delta \bar{z} = \operatorname{Cov}_g\!\left(\bar{w}_g, \bar{z}_g\right) + \mathbb{E}_g\!\left(\bar{w}_g \Delta \bar{z}_g\right)\]

The first term is between-group selection: groups with more cooperators do better as groups. The second is within-group selection: inside any one group, cooperators do worse than the defectors they subsidise. The two terms have opposite signs for a cooperative trait, and the outcome is whichever wins.

This post builds a population where the within-group term is negative in every single group and the population-wide cooperator frequency rises anyway. That is Simpson’s paradox with an evolutionary payload, and it is the honest core of what multilevel selection means.

A group-structured population

Each group has m members. A cooperator pays a cost c and delivers a benefit b split equally among the other m - 1 members of its group, so it receives nothing from its own act. With k cooperators in a group, a cooperator collects benefits from the other k - 1 cooperators and a defector collects from all k.

pcov <- function(x, y) mean(x * y) - mean(x) * mean(y)   # population covariance

make_pop <- function(kvec, m, b, cost, w0 = 1) {
  do.call(rbind, lapply(seq_along(kvec), function(g) {
    k <- kvec[g]
    z   <- c(rep(1, k), rep(0, m - k))                   # 1 = cooperator
    pay <- ifelse(z == 1, -cost + b * (k - 1) / (m - 1), b * k / (m - 1))
    data.frame(group = g, z = z, w = w0 + pay)
  }))
}

pop <- make_pop(kvec = c(2, 8), m = 10, b = 5, cost = 1)
print(aggregate(cbind(z, w) ~ group, pop, mean))
  group   z   w
1     1 0.2 1.8
2     2 0.8 4.2

Two groups of ten, one with two cooperators and one with eight. The cooperator-rich group has a mean fitness of 4.2 against 1.8, which is where the between-group term will come from. The group mean payoff works out to exactly (b - c) times the cooperator fraction, so a group’s productivity is linear in how cooperative it is.

Inside each group the picture is the reverse.

for (g in 1:2) {
  d  <- subset(pop, group == g)
  wc <- mean(d$w[d$z == 1]); wd <- mean(d$w[d$z == 0])
  cat(sprintf("group %d: x = %.2f   cooperator %.4f   defector %.4f   difference %+.6f\n",
              g, mean(d$z), wc, wd, wc - wd))
}
group 1: x = 0.20   cooperator 0.5556   defector 2.1111   difference -1.555556
group 2: x = 0.80   cooperator 3.8889   defector 5.4444   difference -1.555556
cat("algebra says the difference is -c - b/(m-1) =", -1 - 5 / 9, "\n")
algebra says the difference is -c - b/(m-1) = -1.555556 

Defectors out-earn cooperators by the same amount in both groups, and that amount does not depend on the group’s composition at all. A cooperator forgoes c and hands b/(m-1) to each other member, including the defector it is being compared with, so the gap is c + b/(m-1) regardless of k. Within-group selection against cooperation is constant here, which makes the accounting unusually clean.

The partition

price_parts <- function(pop) {
  gm <- aggregate(cbind(z, w) ~ group, pop, mean)
  list(between = pcov(gm$w, gm$z),
       within  = mean(sapply(split(pop, pop$group), function(d) pcov(d$w, d$z))),
       total   = pcov(pop$w, pop$z),
       wbar    = mean(pop$w),
       zbar    = mean(pop$z))
}

p <- price_parts(pop)
cat(sprintf("between %+.6f   within %+.6f   total %+.6f\n", p$between, p$within, p$total))
between +0.360000   within -0.248889   total +0.111111
cat("identity holds:", all.equal(p$total, p$between + p$within), "\n")
identity holds: TRUE 
cat(sprintf("mean fitness %.4f   change in cooperator frequency %+.6f\n",
            p$wbar, p$total / p$wbar))
mean fitness 3.0000   change in cooperator frequency +0.037037
cat(sprintf("cooperator frequency: %.6f before, %.6f after\n",
            p$zbar, p$zbar + p$total / p$wbar))
cooperator frequency: 0.500000 before, 0.537037 after

Between-group selection contributes 0.36 and within-group selection contributes -0.2489, so the total is positive and the cooperator frequency rises from 0.5 to 0.5370. The identity is checked with all.equal rather than asserted, which matters: an identity that has been rearranged by hand is exactly the kind of thing that survives a plausibility check and fails a numeric one.

The paradox made explicit

Track the cooperator fraction inside each group as well as across the population.

library(ggplot2)
after_group <- sapply(1:2, function(g) {
  d <- subset(pop, group == g)
  mean(d$z) * mean(d$w[d$z == 1]) / mean(d$w)
})
before_group <- sapply(1:2, function(g) mean(subset(pop, group == g)$z))
d <- data.frame(
  level = factor(rep(c("group 1", "group 2", "whole population"), each = 2),
                 levels = c("group 1", "group 2", "whole population")),
  time  = rep(c("before", "after"), 3),
  x     = c(before_group[1], after_group[1], before_group[2], after_group[2],
            p$zbar, p$zbar + p$total / p$wbar))
d$time <- factor(d$time, levels = c("before", "after"))
cat(sprintf("group 1: %.6f -> %.6f\ngroup 2: %.6f -> %.6f\nwhole:   %.6f -> %.6f\n",
            d$x[1], d$x[2], d$x[3], d$x[4], d$x[5], d$x[6]))
group 1: 0.200000 -> 0.061728
group 2: 0.800000 -> 0.740741
whole:   0.500000 -> 0.537037
ggplot(d, aes(time, x, group = level, colour = level)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2.6) +
  scale_colour_manual(values = c(`group 1` = "#b5534e", `group 2` = "#c98a3a",
                                 `whole population` = "#275139")) +
  ylim(0, 1) +
  labs(x = NULL, y = "Fraction cooperating", colour = NULL) +
  theme_minimal(base_size = 12) + theme(legend.position = "top")
Three pairs of connected points: two group lines sloping downwards and a population line sloping upwards.
Figure 1: Cooperator fraction before and after one round of selection. Inside both groups the fraction falls; across the population it rises, because the cooperator-rich group contributes more than twice as many offspring.

Group 1 goes from 0.2 to 0.0617 and group 2 from 0.8 to 0.7407. Both fall. The population goes from 0.5 to 0.5370. Nothing has been smuggled in: group 2 simply produces more than twice the offspring of group 1, so the population-wide average shifts towards group 2’s composition faster than either group’s composition erodes.

This is why “cooperators do worse than defectors” and “cooperation is increasing” can both be true statements about the same population at the same moment. Whether they are depends entirely on how much of the trait variance sits between groups rather than within them.

Where the sign flips

Because both terms have closed forms in this model, the condition can be located exactly instead of scanned. Between-group selection is (b - c) times the variance in group composition, and within-group selection is -(c + b/(m-1)) times the average within-group variance.

gm <- aggregate(cbind(z, w) ~ group, pop, mean)
cat("between equals (b - c) Var(x):",
    all.equal(p$between, (5 - 1) * pcov(gm$z, gm$z)), "\n")
between equals (b - c) Var(x): TRUE 
cat("within equals -(c + b/(m-1)) E[x(1-x)]:",
    all.equal(p$within, -(1 + 5 / 9) * mean(gm$z * (1 - gm$z))), "\n")
within equals -(c + b/(m-1)) E[x(1-x)]: TRUE 

Now put two groups at 0.5 - d and 0.5 + d and increase the spread.

components <- function(d, m = 10, b = 5, cost = 1) {
  x   <- c(0.5 - d, 0.5 + d)
  bet <- (b - cost) * pcov(x, x)
  wit <- -(cost + b / (m - 1)) * mean(x * (1 - x))
  c(between = bet, within = wit, total = bet + wit)
}
ds <- seq(0, 0.5, by = 0.005)
cm <- as.data.frame(t(sapply(ds, components)))
cm$d <- ds
print(round(cm[cm$d %in% c(0, 0.1, 0.25, 0.3, 0.4, 0.5), c("d", "between", "within", "total")], 5),
      row.names = FALSE)
    d between   within    total
 0.00    0.00 -0.38889 -0.38889
 0.10    0.04 -0.37333 -0.33333
 0.25    0.25 -0.29167 -0.04167
 0.30    0.36 -0.24889  0.11111
 0.40    0.64 -0.14000  0.50000
 0.50    1.00  0.00000  1.00000
long <- data.frame(d = rep(ds, 3), value = c(cm$between, cm$within, cm$total),
                   term = rep(c("between groups", "within groups", "total"), each = length(ds)))
ggplot(long, aes(d, value, colour = term)) +
  geom_hline(yintercept = 0, colour = "#8d8d80") +
  geom_line(linewidth = 0.9) +
  scale_colour_manual(values = c(`between groups` = "#275139",
                                 `within groups` = "#b5534e", total = "#2f4858")) +
  labs(x = "Between-group spread (d)", y = "Contribution to covariance", colour = NULL) +
  theme_minimal(base_size = 12) + theme(legend.position = "top")
A rising curve for the between-group term, a shallow rising negative curve for the within-group term, and their sum crossing zero at about 0.265.
Figure 2: The two components of change against between-group spread, for two groups of ten with b = 5 and c = 1. Between-group selection grows as the square of the spread while within-group selection weakens, and the total crosses zero near a spread of 0.265.

The crossing point has a closed form. Setting the two terms equal gives a required spread of half the square root of q / ((b - c) + q), where q is the within-group penalty c + b/(m-1).

q     <- 1 + 5 / 9
root  <- uniroot(function(d) components(d)["total"], c(0.01, 0.49), tol = 1e-12)$root
exact <- 0.5 * sqrt(q / (4 + q))
cat(sprintf("numeric root %.8f   closed form %.8f\n", root, exact))
numeric root 0.26457513   closed form 0.26457513
cat("agree:", all.equal(root, exact, tolerance = 1e-6), "\n")
agree: TRUE 
for (m in c(2, 4, 10, 25, 100, 1000)) {
  qm <- 1 + 5 / (m - 1)
  cat(sprintf("m = %4d   within-group penalty %.4f   spread needed %.6f\n",
              m, qm, 0.5 * sqrt(qm / (4 + qm))))
}
m =    2   within-group penalty 6.0000   spread needed 0.387298
m =    4   within-group penalty 2.6667   spread needed 0.316228
m =   10   within-group penalty 1.5556   spread needed 0.264575
m =   25   within-group penalty 1.2083   spread needed 0.240832
m =  100   within-group penalty 1.0505   spread needed 0.228035
m = 1000   within-group penalty 1.0050   spread needed 0.224054
cat("limit for very large groups:", format(0.5 * sqrt(1 / 5), digits = 8), "\n")
limit for very large groups: 0.2236068 

Small groups are harder, not easier, in this parameterisation. In a group of two the cooperator’s entire benefit lands on the single partner it is competing with, so the within-group penalty is c + b, and the required spread is 0.387. In very large groups the benefit is diluted across many members, the penalty falls to c alone, and the required spread drops towards 0.2236. Both limits fall out of the same expression, and neither had to be guessed.

Group size and group variance are not independent in any real system, which is the part this calculation cannot supply. Small groups formed by siblings will have high between-group variance for genetic reasons, and that is where the group and kin accounts meet.

Does it work with many groups?

The identity is not special to two groups. Draw forty groups with random compositions and check both the partition and the closed form.

set.seed(11)
kv   <- sample(0:10, 40, replace = TRUE)
pop2 <- make_pop(kv, m = 10, b = 5, cost = 1)
p2   <- price_parts(pop2)
gm2  <- aggregate(cbind(z, w) ~ group, pop2, mean)
formula_total <- 4 * pcov(gm2$z, gm2$z) - (1 + 5 / 9) * mean(gm2$z * (1 - gm2$z))
cat(sprintf("total from data %.8f   from the formula %.8f\n", p2$total, formula_total))
total from data 0.15187778   from the formula 0.15187778
cat("agree:", all.equal(p2$total, formula_total), "\n")
agree: TRUE 
cat(sprintf("variance in group composition %.6f   change in frequency %+.6f\n",
            pcov(gm2$z, gm2$z), p2$total / p2$wbar))
variance in group composition 0.097275   change in frequency +0.051659

Randomly assembled groups produce enough variance in composition here to push the total positive. That is a warning as much as a result: with b/c = 5 and groups of ten, sampling noise alone can generate the variance that between-group selection needs, and any claim that a real system has group selection must show the variance is larger than random assortment would give.

An honest limit

The partition is an accounting identity, not a mechanism. Every group-structured population can be described this way, whether or not anything biologically interesting is happening, and the same change can be re-described with a different partition and different-looking terms. Choosing the group boundaries chooses the answer. That is not a flaw in the algebra; it is a reason the algebra cannot settle arguments on its own.

The same change also has a kin selection description. Variance in group composition is relatedness by another name, and the between-group term maps onto rb while the within-group term maps onto -c. When people describe multilevel selection and kin selection as rival theories they usually mean two accounting schemes for one process, and the empirical questions (how much variance, at what scale, with what benefit structure) are identical under both.

Finally, this model does one round of selection on fixed groups. Nothing here forms groups, disperses offspring, or lets group composition regenerate, and those steps are where real models earn or lose their conclusions. A population that reshuffles into random groups every generation destroys between-group variance as fast as selection builds it.

References

Price GR 1970. Nature 227(5257):520-521 (10.1038/227520a0)

Wilson DS 1975. Proceedings of the National Academy of Sciences 72(1):143-146 (10.1073/pnas.72.1.143)

Traulsen A, Nowak MA 2006. Proceedings of the National Academy of Sciences 103(29):10952-10955 (10.1073/pnas.0602530103)

Nowak MA 2006. Science 314(5805):1560-1563 (10.1126/science.1133755)

Newsletter

Get new tutorials by email

New R and QGIS tutorials for ecologists, straight to your inbox. No spam; unsubscribe anytime.

By subscribing you agree to receive these emails and confirm your address once. See the privacy policy.