quick_model
In IBMPopSim: Individual Based Model Population Simulation

We take here an initial population, stored in a data frame, composed of $100\,000$ individuals marked by their gender (encoded by a Boolean characteristic):

library(IBMPopSim)

pop <- population(EW_pop_14$sample)

The second step is to define the model parameters' list:

params <- list("alpha" = 0.008, "beta" = 0.02, "p_male" = 0.51,
               "birth_rate" = stepfun(c(15,40), c(0,0.05,0)))

The last step is to defined the events that can occur in the population, here birth and death events:

death_event <- mk_event_individual(type = "death",
                  intensity_code = "result = alpha * exp(beta * age(I, t));")

birth_event <- mk_event_individual(type = "birth", 
                  intensity_code = "result = birth_rate(age(I,t));",
                  kernel_code = "newI.male = CUnif(0, 1) < p_male;")

Note that these events contain some C++ statements that depend (implicitly) on the previously declared parameters in variable params.

Finally, the model is created by calling the function mk_model. A C++ source code is obtained from the events and parameters, then compiled using the sourceCpp function of the Rcpp package.

model <- mk_model(characteristics = get_characteristics(pop),
                  events = list(death_event, birth_event),
                  parameters = params)

In order to simulate a random trajectory of the population until a given time T bounds on the events intensities have to be specified:

a_max <- 115
events_bounds = c("death" = params$alpha * exp(params$beta * a_max),
                  "birth" = max(params$birth_rate))

Then, the function popsim can be called:

sim_out <- popsim(model, pop, events_bounds, params,
                  age_max = a_max, time = 30)

The data frame sim_out$population contains the information (date of birth, date of death, gender) on individuals who lived in the population over the period $[0,30]$. Functions of the package allows to provide aggregated information on the population.

pop_out <- sim_out$population
head(pop_out)
female_pop <- pop_out[pop_out$male==FALSE, ]
age_pyramid(female_pop, ages = 85:90, time = 30)
Dxt <- death_table(female_pop, ages = 85:90, period = 20:30)
Ext <- exposure_table(female_pop, ages = 85:90, period = 20:30)

Note that parameters of the model can be changed without recompiling the model.

params$beta <- 0.01

# Update death event bound:
events_bounds["death"] <- params$alpha * exp(params$beta * a_max)

sim_out <- popsim(model, pop, events_bounds, params,
                  age_max = a_max, time = 30)