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

\[ p(t)=\lambda e^{-\lambda t} \]

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:

downloading and extracting the examples archive from the neworder releases page, or
pulling the docker image virgesmith/neworder:latest (more here

Input¶

Here's the code that sets up and runs the model:

./examples/chapter1/model.py

"""
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
from person import People  # ty:ignore[unresolved-import]

import neworder

# 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(f"Age {age:.0f} survival rate = {model.alive(age):.1%}")

Implementation¶

Each individual is an instance of the Person class:

./examples/chapter1/person.py

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

./examples/chapter1/person.py

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(f"created {n} individuals")

The single timestep records each person's time of death

./examples/chapter1/person.py

    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:

./examples/chapter1/person.py

    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(f"Life expectancy = {sample_le:.2f} years (sampling error={error:.2f} years)")

Then, from the model script, this function is called, displaying the proportion of the cohort that are still alive at 10-year intervals:

./examples/chapter1/person.py

    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%

and clearly a constant mortality rate isn't realistic as we see far more deaths at younger ages, and far less at older ages, than would be expected. The example mortality introduces a model with a time-dependent mortality hazard rate and shows how the framework can very efficiently model this.