Chapter 1¶
This example is based on the example in the introductory chapter of Microsimulation and Population Dynamics. This is a simple example - a basic cohort model using continuous-time case-based simulation of mortality (only) for a homogeneous population, using a constant mortality hazard rate. In other words, age at death is sampled from an exponential distribution
which has a mean, i.e. life expectancy, of \(\mu=1/\lambda\).
The neworder implementation is as direct a port of the MODGEN model, as far as possible.
Examples Source Code
The code for all the examples can be obtained by either:
Input¶
Here's the code that sets up and runs the model:
"""
Chapter 1
This is a direct neworder cover version of the Basic Cohort Model from the Belanger & Sabourin book
See https://www.microsimulationandpopulationdynamics.com/
"""
import numpy as np
import neworder
from person import People
# neworder.verbose() # uncomment for detailed output
# "An arbitrarily selected value, chosen to produce a life expectancy of about 70 years."
# (see the "mortality" example for a more realistic model)
mortality_hazard = 0.014
population_size = 100000
model = People(mortality_hazard, population_size)
neworder.run(model)
# now we can sample the population generated by the model to see the proportion of deaths at (arbitrarily) 10 year intervals
for age in np.linspace(10.0, 100.0, 10):
neworder.log("Age %.0f survival rate = %.1f%%" % (age, model.alive(age) * 100.0))
Implementation¶
Each individual is an instance of the Person class:
class Person():
"""
MODGEN equivalent: actor Person {...}
Represents a single individual
"""
def __init__(self, mortality_hazard: float) -> None:
""" MODGEN equivalent: Person::Start() """
self.alive = True
# MODGEN would automatically create time and age, though they are not needed for this simple example
self.mortality_hazard = mortality_hazard
self.time_mortality = neworder.time.NEVER # to be computed later
def finish(self) -> None:
""" MODGEN equivalent: Person::Finish() """
# nothing required here
def state(self, t: float) -> bool:
""" Returns the person's state (alive/dead) at age t """
return self.time_mortality > t
def time_mortality_event(self, mc: neworder.MonteCarlo) -> float:
""" MODGEN equivalent: TIME Person::timeMortalityEvent() """
if neworder.time.isnever(self.time_mortality):
self.time_mortality = mc.stopping(self.mortality_hazard, 1)[0]
return self.time_mortality
def mortality_event(self) -> None:
""" MODGEN equivalent: void Person::MortalityEvent() """
# NB this is not used in this implementation
self.alive = False
And the People model contains an array of Persons
class People(neworder.Model):
""" A model containing an aggregration of Person objects """
def __init__(self, mortality_hazard: float, n: int) -> None:
# initialise base model with a nondeterministic seed results will vary (slightly)
super().__init__(neworder.NoTimeline(), neworder.MonteCarlo.nondeterministic_stream)
# initialise population
self.population = [Person(mortality_hazard) for _ in range(n)]
neworder.log("created %d individuals" % n)
The single timestep records each person's time of death
def step(self) -> None:
# sample each person's age at death.
# (this is not an efficient implementation when everyone has the same hazard rate)
[p.time_mortality_event(self.mc) for p in self.population]
And when the timeline has reached the end, the finalise method compares the mean of the sampled times of death with the expected value:
def finalise(self) -> None:
# compute mean sampled life expectancy against theoretical
sample_le = sum([p.time_mortality for p in self.population]) / len(self.population)
actual_le = 1.0 / self.population[0].mortality_hazard
error = sample_le - actual_le
neworder.log("Life expectancy = %.2f years (sampling error=%.2f years)" % (sample_le, error))
Then, from the model script, this function is called, displaying the proportion of the cohort that are still alive at 10-year intervals:
def alive(self, t: float) -> float:
return np.mean([p.state(t) for p in self.population])
Output¶
Some example output:
[py 0/1] created 100000 individuals
[py 0/1] Life expectancy = 71.60 years (sampling error=0.17 years)
[py 0/1] Age 10 survival rate = 86.9%
[py 0/1] Age 20 survival rate = 75.6%
[py 0/1] Age 30 survival rate = 65.8%
[py 0/1] Age 40 survival rate = 57.3%
[py 0/1] Age 50 survival rate = 49.8%
[py 0/1] Age 60 survival rate = 43.3%
[py 0/1] Age 70 survival rate = 37.5%
[py 0/1] Age 80 survival rate = 32.7%
[py 0/1] Age 90 survival rate = 28.3%
[py 0/1] Age 100 survival rate = 24.7%