In nick-gauthier/Silvanus: Multi-Level Simulation of Roman Land Use

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 6,
  fig.asp = 0.618
)

library(dplyr)
library(purrr)
library(tidyr)
library(ggplot2)
library(patchwork)
library(Silvanus)
library(furrr)

Introduction

Get the life table here and talk about the basic

First we'll create a test population of 15 25-year-olds.

pop1 <- create_individuals(15)

Unbounded growth

Then run a 500 year simulation from that initial population.

set.seed(1000) # set seed for reproducibility
nsim <- 500

demographic_sim <- tibble(time = 0:nsim,
                          age = accumulate(1:nsim, ~population_dynamics(.x), .init = pop1),
                          population = map_dbl(age, length)) %>%
  mutate(dependency_ratio = map_dbl(age,  ~sum(!between(.x, 15, 65)) / sum(between(.x, 15, 65))) * 100)

Look at the results.

Dependency ratio

a <-   ggplot(demographic_sim, aes(time, dependency_ratio))+
  geom_line() +
  labs(x = 'Year', y = 'Dependency ratio') +
  theme_bw()
d1 <- demographic_sim %>%
  filter(time %in% c(0, 100, 200, 300, 400, 500)) %>%
  unnest(cols = age) %>%
  mutate(time = paste('Year', time)) %>%
  group_by(time) %>%
  mutate(median_age = median(age))
b <- ggplot(d1, aes(age)) +
  geom_histogram(binwidth = 1, center = 0) +
  geom_vline(aes(xintercept = median_age), color = 'red', linetype = 2) + # find a way to get median age to vary with facet
  facet_wrap(~time, scales = 'free_y') +
  labs(x = 'Age', y = 'Count') +
 # labs(title = paste('Evolution of the age distribution over', nsim, 'years'), subtitle = 'Century snapshots') +
  theme_bw()

a / b + plot_annotation(tag_levels = 'A')

Growth rates

plan(multisession)

r_test <- future_map(1:100, ~ 
                       tibble(time = 0:600,
                          age = accumulate(1:600, ~population_dynamics(.x), .init = pop1),
                          population = map_dbl(age, length)) %>%
                       select(-age), .progress = TRUE) %>%
  bind_rows(.id = 'rep')

Fit a poisson GLM with a log link function and a varying intercept to represent the different effective initial conditions of each simulation run.

m1 <- glm(population ~ time + rep, data = r_test, family = poisson(link = 'log'))

From the GLM, we see that the intrinsic growth rate is r round(coefficients(m1)['time'] * 100, 2)% per year.

ggplot(r_test, aes(time, population)) +
  geom_line(aes( group = rep), alpha = .1) +
  geom_smooth(method = 'lm', color = 'red') +
  scale_y_log10() +
  theme_classic()

Food-limited growth

Now we can look at a food-limited population.

plan(multisession)
demographic_sim2  <- expand_grid(rep = 1:25, 
                                 food_ratio = seq(0.7, 1, by = 0.05)) %>%
  mutate(sim = future_map(food_ratio, function(x) accumulate(1:nsim, ~population_dynamics(.x, food_ratio_c = x), .init = pop1), .progress = TRUE),
         sim = map(sim, ~tibble(time = 0:nsim,
                          age = .,
                          population = map_dbl(age, length)) %>%
                     mutate(dependency_ratio = map_dbl(age,  ~sum(!between(.x, 15, 65))/sum(between(.x, 15, 65))) * 100)))

demographic_sim2 %>% unnest %>%  filter(time == nsim) %>%
  ggplot(aes(food_ratio, population, group = food_ratio)) +
  geom_boxplot() +
  theme_classic()

demographic_sim2 %>% unnest %>%  filter(time == nsim) %>%
  ggplot(aes(food_ratio, dependency_ratio, group = food_ratio)) +
  geom_boxplot() +
  theme_classic()

demographic_sim2 %>%
  unnest %>%
 # filter(food_ratio >= .8) %>%
  ggplot(aes(time, population, group =interaction(rep, food_ratio), color = food_ratio)) +
  geom_line(alpha = .5) +
  #scale_y_continuous(trans = 'log') +
  scale_color_viridis_c() +
  theme_bw()

demographic_sim2 %>%
  unnest %>%
 filter(food_ratio >= .8) %>%
  ggplot(aes(time, population, group =interaction(rep, food_ratio), color = food_ratio)) +
  geom_line(alpha = .5) +
  scale_y_continuous(trans = 'log') +
  scale_color_viridis_c() +
  theme_bw()

Drought response

nsim <- 500
rn <- c(rep(1, 300), rep(0.5, 30), rep(1, 170))
demographic_sim3 <- tibble(time = 0:nsim,
                          age = accumulate(rn, ~population_dynamics(.x, food_ratio_c = .y), .init = pop1),
                          population = map_dbl(age, length)) %>%
  mutate(dependency_ratio = map_dbl(age,  ~sum(!between(.x, 15, 65))/sum(between(.x, 15, 65))) * 100)

demographic_sim3 <- run_sim(t1, food_ratio = rn, nsim = nsim)

ggplot(demographic_sim3, aes(time, population)) + 
  geom_ribbon(aes(xmin = 300, xmax = 330), alpha = 0.2) +
  geom_line() +
  theme_classic()


ggplot(demographic_sim3, aes(time, dependency_ratio)) + 
  geom_ribbon(aes(xmin = 300, xmax = 330), alpha = 0.2) +
  geom_line() +
  theme_classic()

create_individuals <- function(n) {
  tibble(age = 0:119) %>%
   mutate(n = if_else(age != 15, 0L, as.integer(n)),
         food_ratio = 1) %>%
       left_join(life_table, by = 'age')
}
t1 <- create_individuals(25)

population_dynamics2 <- function(individuals, food_ratio_c) {
  individuals %>%
    mutate(food_ratio = food_ratio_c) %>%
    reproduce2() %>%
    die2() %>%
    mutate(n = lag(n, default = 0))# happy birthday!
}

reproduce2 <- function(individuals) {
    individuals %>%
    mutate(fertility_reduction = pgamma(food_ratio, shape = fertility_shape, scale = fertility_scale),
         births = rbinom(n(), n, fertility_rate / 2 * fertility_reduction),
         n = if_else(age != 0, n, n + sum(births))) %>%
    select(-fertility_reduction, -births)
}

die2 <- function(individuals) {
  individuals %>%
        mutate(survival_reduction = pgamma(pmin(1, food_ratio), shape = survival_shape, scale = survival_scale),
           survivors = rbinom(n(), n, survival_rate * survival_reduction),
           n = survivors) %>%
    select(-survival_reduction, -survivors, -food_ratio)
}

postprocess <- function(individuals) {
  individuals %>%
  bind_rows(.id = 'time') %>%
  mutate(time = as.numeric(time),
         worker = between(age, 15, 65)) %>%
  group_by(time) %>%
  count(worker, wt = n) %>%
  ungroup %>%
  spread(worker, n) %>%
  mutate(population = `FALSE` + `TRUE`,
         dependency_ratio = `FALSE` / `TRUE` * 100) %>%
    select(-`FALSE`, -`TRUE`)
}

run_sim <- function(individuals, nsim, food_ratio) {
  #length(food_ratio) should be either nsim or 1
  if(length(food_ratio) == 1) food_ratio <- rep(food_ratio, nsim)

  accumulate(food_ratio, ~population_dynamics2(.x, food_ratio_c = .y), .init = individuals) %>%
  postprocess()
}

demsim <- run_sim(t1, nsim = 500, food_ratio = 1)



ggplot(demsim, aes(time, population)) + geom_line() +
  scale_y_log10()

ggplot(demsim, aes(time, dependency_ratio)) + geom_line()

r_test <- future_map_dfr(1:100, 
                     ~accumulate(1:500, ~population_dynamics2(., food_ratio = 1), .init = t1) %>%
                       postprocess(), 
                     .id = 'rep', .progress = TRUE)