repet time to extinction experiment

Chapter_8.Rmd

@@ -135,3 +135,322 @@ par(old_par)
 ```
 
 If the  skill level is transmitted culturally  from these skilled individual to the next generation, we can then add different type of social transmission process to then get prediction about on the spread (and retention) of skills within the population.
+
+## Random copying and the loss of highly skilled people
+
+We now draw one initial skill for each person from the gamma distribution above,
+keeping the same populations of 500 and 1,000 people and the high-skill threshold
+of 6. These skills form generation 0 and are generated only once.
+
+In each subsequent generation, every person copies the skill of a randomly
+chosen person in the immediately preceding generation. We sample **with
+replacement**, so the same person can be copied several times, while others may
+not be copied at all. Each population keeps its original size throughout the
+simulation. Copying is exact, and everyone has the same chance of being chosen,
+regardless of their skill.
+
+We continue sampling until both populations have no one above the threshold.
+Random copying can also leave everyone in a population above the threshold;
+that state is permanent under exact copying, so we treat it as an alternative
+outcome. The simulation stops once both populations have reached one of these
+outcomes, with a limit of 10,000 generations if either population remains mixed.
+
+```{r skill-random-copying}
+set.seed(42)
+max_sampling_generations <- 10000
+sampling_generations <- 0L
+
+initial_skill_populations <- lapply(population_sizes, function(size) {
+  rgamma(size, shape = skill_shape, scale = skill_scale)
+})
+skill_populations <- initial_skill_populations
+
+# Include generation 0 so we can compare the initial and inherited skills.
+skilled_counts <- matrix(
+  0L,
+  nrow = max_sampling_generations + 1,
+  ncol = length(population_sizes),
+  dimnames = list(0:max_sampling_generations, paste("N =", population_sizes))
+)
+current_skilled_counts <- vapply(
+  skill_populations,
+  function(skills) sum(skills > high_skill_cutoff),
+  integer(1)
+)
+skilled_counts[1, ] <- current_skilled_counts
+
+# A mixed population can still gain or lose people above the threshold.
+while (
+  sampling_generations < max_sampling_generations &&
+  any(current_skilled_counts > 0 & current_skilled_counts < population_sizes)
+) {
+  # Each new population copies from its own previous generation.
+  skill_populations <- lapply(skill_populations, function(skills) {
+    sample(skills, size = length(skills), replace = TRUE)
+  })
+  sampling_generations <- sampling_generations + 1L
+  current_skilled_counts <- vapply(
+    skill_populations,
+    function(skills) sum(skills > high_skill_cutoff),
+    integer(1)
+  )
+  skilled_counts[sampling_generations + 1, ] <- current_skilled_counts
+}
+
+# Discard unused rows after the stopping condition is reached.
+skilled_counts <- skilled_counts[seq_len(sampling_generations + 1), , drop = FALSE]
+```
+
+The histograms show the skills remaining after `r sampling_generations`
+generations. Both panels use the same bins and axes. The dashed line marks the
+threshold, and the darker bars count people whose skills are strictly above it.
+
+```{r skill-random-copying-distribution, fig.width=9, fig.height=4.5, fig.cap=paste("Skill distributions after", sampling_generations, "generations of random copying with replacement. Population sizes stay constant, but the frequencies of inherited skills change.")}
+# Make the threshold a bin boundary so the shaded counts match the cutoff.
+skill_breaks <- sort(unique(c(
+  seq(0, ceiling(max(unlist(initial_skill_populations))), by = 0.5),
+  high_skill_cutoff
+)))
+final_skill_histograms <- lapply(skill_populations, function(skills) {
+  hist(skills, breaks = skill_breaks, plot = FALSE)
+})
+skill_histogram_y_max <- 1.15 * max(vapply(
+  final_skill_histograms,
+  function(skill_histogram) max(skill_histogram$counts),
+  numeric(1)
+))
+
+old_par <- par(mfrow = c(1, 2), mar = c(4.5, 4.5, 3.5, 1))
+
+for (i in seq_along(population_sizes)) {
+  plot(
+    final_skill_histograms[[i]],
+    freq = TRUE,
+    col = ifelse(
+      final_skill_histograms[[i]]$mids > high_skill_cutoff,
+      tail_cols[i],
+      body_cols[i]
+    ),
+    border = "white",
+    xlim = range(skill_breaks),
+    ylim = c(0, skill_histogram_y_max),
+    main = paste("Population:", format(population_sizes[i], big.mark = ",")),
+    xlab = "Skill level",
+    ylab = "Number of people",
+    las = 1
+  )
+  abline(v = high_skill_cutoff, lty = 3, col = tail_cols[i])
+  mtext(
+    paste(
+      skilled_counts[sampling_generations + 1, i],
+      "people with skill >", high_skill_cutoff
+    ),
+    side = 3,
+    line = 0.3,
+    cex = 0.85,
+    col = line_cols[i]
+  )
+}
+
+par(old_par)
+```
+
+We can also follow the number of people above the same threshold at each
+generation. Unlike the expected counts in the earlier example, these are counts
+of individuals in the simulated populations.
+
+```{r skill-random-copying-counts, fig.width=7.5, fig.height=4.5, fig.cap=paste("The number of people with skill above", high_skill_cutoff, "across", sampling_generations, "generations. Generation 0 is the initial population; subsequent changes result only from random copying.")}
+matplot(
+  0:sampling_generations,
+  skilled_counts,
+  type = "l",
+  lty = c(1, 2),
+  lwd = 2,
+  col = line_cols,
+  ylim = c(0, max(1, skilled_counts) * 1.2),
+  xlab = "Generation",
+  ylab = paste("People with skill >", high_skill_cutoff),
+  xaxt = "n",
+  las = 1
+)
+axis(1, at = pretty(c(0, sampling_generations)))
+legend(
+  "topright",
+  legend = paste("Population:", format(population_sizes, big.mark = ",")),
+  col = line_cols,
+  lty = c(1, 2),
+  lwd = 2,
+  bty = "n"
+)
+```
+
+```{r skill-random-copying-summary}
+skill_outcome_generations <- vapply(seq_along(population_sizes), function(i) {
+  reached_outcome <- which(
+    skilled_counts[, i] == 0 | skilled_counts[, i] == population_sizes[i]
+  )
+  if (length(reached_outcome) > 0) {
+    reached_outcome[1] - 1L
+  } else {
+    sampling_generations
+  }
+}, integer(1))
+skill_outcomes <- ifelse(
+  current_skilled_counts == 0,
+  "No one above threshold",
+  ifelse(
+    current_skilled_counts == population_sizes,
+    "Everyone above threshold",
+    "Generation limit reached"
+  )
+)
+
+knitr::kable(
+  data.frame(
+    Population = population_sizes,
+    `Initial count` = skilled_counts[1, ],
+    `Final count` = current_skilled_counts,
+    Generation = skill_outcome_generations,
+    Outcome = skill_outcomes,
+    row.names = NULL,
+    check.names = FALSE
+  ),
+  caption = paste(
+    "Number of people with skill above", high_skill_cutoff,
+    "and the generation when each outcome was reached."
+  )
+)
+```
+
+The number of highly skilled people can rise or fall purely by chance. This is
+**cultural drift**: random copying changes the frequencies of skills without a
+preference for higher or lower skill. In expectation, the number above the
+threshold in the next generation equals the current number, but an individual
+run can move away from its starting point. Skills that disappear cannot return
+in this model because there is no innovation or copying error.
+
+## So what happen _in average_?
+
+If we repeat the experiment 200 times for each population size. we can record when the number of skilled worker become null. 
+
+
+```{r skill-extinction-repeats}
+set.seed(2026)
+extinction_repeats <- 200
+
+extinction_results <- expand.grid(
+  population_size = population_sizes,
+  run = seq_len(extinction_repeats)
+)
+extinction_results$initial_skilled <- NA_integer_
+extinction_results$final_skilled <- NA_integer_
+extinction_results$generations <- NA_integer_
+extinction_results$extinction_time <- NA_integer_
+extinction_results$outcome <- NA_character_
+
+for (r in seq_len(nrow(extinction_results))) {
+  repeat_size <- extinction_results$population_size[r]
+  repeat_skills <- rgamma(repeat_size, shape = skill_shape, scale = skill_scale)
+  repeat_skilled_count <- sum(repeat_skills > high_skill_cutoff)
+  extinction_results$initial_skilled[r] <- repeat_skilled_count
+  repeat_generation <- 0L
+
+  while (
+    repeat_generation < max_sampling_generations &&
+    repeat_skilled_count > 0 && repeat_skilled_count < repeat_size
+  ) {
+    repeat_skills <- sample(repeat_skills, size = repeat_size, replace = TRUE)
+    repeat_generation <- repeat_generation + 1L
+    repeat_skilled_count <- sum(repeat_skills > high_skill_cutoff)
+  }
+
+  extinction_results$final_skilled[r] <- repeat_skilled_count
+  extinction_results$generations[r] <- repeat_generation
+  if (repeat_skilled_count == 0) {
+    extinction_results$extinction_time[r] <- repeat_generation
+    extinction_results$outcome[r] <- "Extinction"
+  } else if (repeat_skilled_count == repeat_size) {
+    extinction_results$outcome[r] <- "Everyone above threshold"
+  } else {
+    extinction_results$outcome[r] <- "Generation limit reached"
+  }
+}
+
+# Only observed extinctions have an extinction time; other outcomes remain NA.
+extinction_times <- lapply(population_sizes, function(size) {
+  times <- extinction_results$extinction_time[
+    extinction_results$population_size == size
+  ]
+  times[!is.na(times)]
+})
+
+extinction_summary <- data.frame(
+  Population = population_sizes,
+  Runs = extinction_repeats,
+  Extinctions = lengths(extinction_times),
+  `All above threshold` = vapply(population_sizes, function(size) {
+    sum(extinction_results$population_size == size &
+          extinction_results$outcome == "Everyone above threshold")
+  }, integer(1)),
+  `At generation limit` = vapply(population_sizes, function(size) {
+    sum(extinction_results$population_size == size &
+          extinction_results$outcome == "Generation limit reached")
+  }, integer(1)),
+  `Median extinction time` = vapply(extinction_times, function(times) {
+    if (length(times) > 0) median(times) else NA_real_
+  }, numeric(1)),
+  check.names = FALSE
+)
+```
+
+```{r skill-extinction-boxplot, fig.width=3.5, fig.height=7.5, fig.cap=paste("Time to extinction among runs with observed extinction, from", extinction_repeats, "runs per population size. ")}
+old_par <- par(mar = c(6, 4.5, 1, 1))
+
+if (any(lengths(extinction_times) > 0)) {
+  boxplot(
+    extinction_times,
+    names = format(population_sizes, big.mark = ",", trim = TRUE),
+    col = body_cols,
+    border = line_cols,
+    outpch = 16,
+    outcex = 0.7,
+    ylim = range(0, unlist(extinction_times)),
+    range=0,
+    xlab = "",
+    ylab = "Generations until extinction",
+    lig="y",
+    las = 1
+  )
+  mtext(
+    paste0(lengths(extinction_times), "/", extinction_repeats, " extinctions"),
+    side = 1,
+    at = seq_along(population_sizes),
+    line = 2.3,
+    cex = 0.9
+  )
+  mtext("Population size", side = 1, line = 4.3)
+} else {
+  plot.new()
+  text(0.5, 0.5, "No extinctions were observed within the generation limit.")
+}
+
+par(old_par)
+```
+
+```{r skill-extinction-summary}
+knitr::kable(
+  extinction_summary,
+  digits = 1,
+  caption = paste(
+    "Outcomes across", extinction_repeats,
+    "runs per population size. Median extinction times are in generations",
+    "and include only observed extinctions."
+  )
+)
+```
+
+The spread of extinction times shows why one run is insufficient to characterise
+the loss of highly skilled people. Larger populations can retain these skills
+for longer, but outcomes still vary because copying is random. These boxplots
+describe the timing of loss **conditional on extinction being observed**; the
+table also shows how often extinction occurred among all the runs.