431 lines
14 KiB
Text
431 lines
14 KiB
Text
# Chapter 8
|
|
|
|
## Population size and highly skilled individuals
|
|
|
|
We can use some theoretical distribution to represent the skills of people in a population
|
|
Using a gamma distribution here we create a scenario where most people are concentrated at lower skill levels, while a longer tail extends towards high skill.
|
|
The high-skill threshold is fixed, so increasing population size.
|
|
If people copu only the 2% most skliied individuals, then a bigger population wil mathematically increases the number of skilled individual
|
|
```{r skill-population-calculation}
|
|
population_sizes <- c(500, 1000)
|
|
skill_shape <- 2
|
|
skill_scale <- 1.5
|
|
high_skill_cutoff <- 6
|
|
|
|
high_skill_proportion <- pgamma(
|
|
high_skill_cutoff,
|
|
shape = skill_shape,
|
|
scale = skill_scale,
|
|
lower.tail = FALSE
|
|
)
|
|
|
|
high_skill_counts <- round(population_sizes * high_skill_proportion)
|
|
|
|
data.frame(
|
|
Population = population_sizes,
|
|
`High-skill proportion` = round(high_skill_proportion, 3),
|
|
`Expected high-skill individuals` = high_skill_counts,
|
|
check.names = FALSE
|
|
)
|
|
```
|
|
|
|
```{r skill-population-figure, fig.width=9, fig.height=4.5, fig.cap="The same skill distribution in two populations. The blue tail marks individuals above the high-skill threshold."}
|
|
skill_x <- seq(
|
|
0,
|
|
qgamma(0.999, shape = skill_shape, scale = skill_scale),
|
|
length.out = 600
|
|
)
|
|
tail_x <- skill_x[skill_x >= high_skill_cutoff]
|
|
|
|
skill_density <- dgamma(
|
|
skill_x,
|
|
shape = skill_shape,
|
|
scale = skill_scale
|
|
)
|
|
tail_density <- dgamma(
|
|
tail_x,
|
|
shape = skill_shape,
|
|
scale = skill_scale
|
|
)
|
|
common_y_max <- max(population_sizes) * max(skill_density) * 1.18
|
|
|
|
line_cols <- c("#244A68", "#355E3B")
|
|
body_cols <- c("#DCE8F0", "#D4DFC8")
|
|
tail_cols <- c("#1769AA", "#3E713F")
|
|
|
|
old_par <- par(
|
|
mfrow = c(1, 2),
|
|
mar = c(4.5, 4.5, 3, 1),
|
|
oma = c(0, 1.5, 0, 0)
|
|
)
|
|
|
|
for (i in seq_along(population_sizes)) {
|
|
people_density <- population_sizes[i] * skill_density
|
|
high_skill_density <- population_sizes[i] * tail_density
|
|
|
|
plot(
|
|
skill_x,
|
|
people_density,
|
|
type = "n",
|
|
xlim = range(skill_x),
|
|
ylim = c(0, common_y_max),
|
|
axes = FALSE,ylab="",xlab = ""
|
|
)
|
|
|
|
polygon(
|
|
c(min(skill_x), skill_x, max(skill_x)),
|
|
c(0, people_density, 0),
|
|
col = body_cols[i],
|
|
border = NA
|
|
)
|
|
polygon(
|
|
c(high_skill_cutoff, tail_x, max(tail_x)),
|
|
c(0, high_skill_density, 0),
|
|
col = tail_cols[i],
|
|
border = NA
|
|
)
|
|
lines(skill_x, people_density, lwd = 2.2, col = line_cols[i])
|
|
abline(v = high_skill_cutoff, lty = 3, col = tail_cols[i])
|
|
|
|
axis(
|
|
1,
|
|
at = c(
|
|
qgamma(0.1, shape = skill_shape, scale = skill_scale),
|
|
qgamma(0.93, shape = skill_shape, scale = skill_scale)
|
|
),
|
|
labels = c("Low skill", "High skill"),
|
|
lwd = 0,
|
|
lwd.ticks = 1
|
|
)
|
|
axis(2, las = 1, lwd = 0, lwd.ticks = 1)
|
|
mtext("Skill level", side = 1, line = 2.5)
|
|
arrows(
|
|
min(skill_x), 0,
|
|
max(skill_x), 0,
|
|
length = 0.08,
|
|
lwd = 1.2,
|
|
xpd = NA
|
|
)
|
|
arrows(
|
|
min(skill_x), 0,
|
|
min(skill_x), common_y_max,
|
|
length = 0.08,
|
|
lwd = 1.2,
|
|
xpd = NA
|
|
)
|
|
|
|
title(main = paste(
|
|
"Population:",
|
|
format(population_sizes[i], big.mark = ",")
|
|
))
|
|
text(
|
|
high_skill_cutoff + 0.38 * (max(skill_x) - high_skill_cutoff),
|
|
people_density[which.max(skill_density)] * 0.62,
|
|
labels = paste(
|
|
format(high_skill_counts[i], big.mark = ","),
|
|
"high-skill\nindividuals"
|
|
),
|
|
col = tail_cols[i],
|
|
font = 2
|
|
)
|
|
}
|
|
|
|
mtext("Number of people", side = 2, outer = TRUE, line = 0.2)
|
|
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)
|
|
```
|
|
|