1 changed files with 0 additions and 294 deletions
--- a/Chapter_8.Rmd
+++ b/Chapter_8.Rmd
@ -135,297 +135,3 @@ 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=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 = c(10,5000),
    range=0,
    xlab = "Population size",
    ylab = "Generations until extinction",
    log="y",
    las = 1
  )
 } else {
  plot.new()
  text(0.5, 0.5, "No extinctions were observed within the generation limit.")
 }
 par(old_par)
 ```