Epidemic forecasting in Python¶

Simulation parameters¶

Default values for the particle filter parameters are provided:

epifx.default_params(p_exp, px_scale, model=epifx.SEIR)¶

The default simulation parameters.

Parameters:	p_exp – The daily probability of seeding an exposure event in a naive population. px_scale – The ratio of particles to seeded exposures. model – The infection model.

Note that epifx uses an independent pseudo-random number generator (PRNG) for sampling from model priors and adding stochastic noise to model parameters and flows (params['epifx']['rnd']) as distinct from the pypfilt PRNG instance (params['rnd']). This PNRG can be initialised with a custom seed:

epifx.init_prng(params, seed)¶

Initialise the epifx PRNG instance (params['epifx']['rnd']).

Parameters:	params – The simulation parameters. seed – The PRNG seed; see the numpy.random.RandomState documentation for details.

The behaviour of epifx 0.1.x (one common PRNG instance) can be recovered by setting the following parameter:

params['epifx']['independent_prng'] = False

Temporal forcing¶

The force of infection can be subject to temporal forcing (e.g., by climatic factors such as absolute humidity, or social factors such as media coverage) as mediated by the model parameter \(\sigma\) (see Infection models). This requires the simulation parameters to include a function that maps datetime.datetime instances to scalar values:

epifx.daily_forcing(filename, date_fmt=u'%Y-%m-%d')¶

Return a temporal forcing look-up function, which should be stored in params['epifx']['forcing'] in order to enable temporal forcing.

Parameters:	filename – A file that contains two columns separated by whitespace, the column first being the date and the second being the force of the temporal modulation. Note that the first line of this file is assumed to contain column headings and will be ignored. date_fmt – The format in which dates are stored.

Infection models¶

class epifx.SEIR¶

An SEIR compartment model for a single circulating influenza strain, under the assumption that recovered individuals are completely protected against reinfection.

\[\begin{split}\frac{dS}{dt} &= - \alpha \cdot S^\eta I \\[0.5em] \frac{dE}{dt} &= \alpha \cdot S^\eta I - \beta \cdot E \\[0.5em] \frac{dI}{dt} &= \beta \cdot E - \gamma I \\[0.5em] \frac{dR}{dt} &= \gamma I\end{split}\]

Parameter	Meaning
\(\alpha\)	Force of infection
\(\beta\)	Incubation period (day -1)
\(\gamma\)	Infectious period (day -1)
\(\eta\)	Inhomogeneous social mixing
\(\sigma\)	Temporal forcing coefficient

The force of infection can be subject to temporal forcing \(F(t)\), as mediated by \(\sigma\):

\[\alpha'(t) = \alpha \cdot (1 + \sigma \cdot F(t))\]

Note that this requires the forcing function \(F(t)\) to be defined in the simulation parameters (see the Temporal Forcing documentation).

static priors(params)¶

Return a dictionary of model parameter priors.

Parameters:	params – Simulation parameters.

Observation models¶

Observation models comprise a log-likelihood function (log_llhd) that relates observations to particles, and an expectation function (expect) that relates particles to observations:

def log_llhd(params, op, obs, pr_indiv, curr, hist):
    """Calculate the log-likelihood of obtaining one or more observations
    from each particle.

    :param params: The simulation parameters.
    :param op: The observation model parameters.
    :param obs: The list of observations for the current time-step.
    :param pr_indiv: A 2 x N matrix (for N particles) of the probability
        of an individual being uninfected and infected, respectively.
    :param curr: The particle state vectors.
    :param hist: The particle state histories, indexed by observation period.
    """
    return ...

def expect(params, op, period, pr_inf, prev, curr):
    """Determine the expected value for a given infection probability.

    :param params: The observation model parameters.
    :param op: The observation model parameters.
    :param period: The duration of the observation period (in days).
    :param pr_inf: The probability of an individual becoming infected
        during the observation period.
    :param prev: The state vectors at the start of the observation period.
    :param curr: The state vectors at the end of the observation period.
    """
    return ...

Observation models are identified by unit, and must be stored in the parameters dictionary as shown below:

params = epifx.default_params(...)
obs_unit = 'Some identifier'
params['obs'][obs_unit] = {
    'log_llhd_fn': log_llhd,
    'expect_fn': expect,
    # Add observation model parameters here.
}

The epifx.obs module provides generic observation models for both count data and fractional data, as well as functions for reading observations from disk.

Forecast summaries¶

epifx.summary.make(params, all_obs, extra_tbls=None, pkgs=None, **kwargs)¶

A convenience function that collects all of the summary statistics defined in the pypfilt.summary and epifx.summary modules.

Parameters:

params – The simulation parameters. all_obs – A list of all observations. extra_tbls – A list of extra summary statistic tables to include. pkgs – A dictionary of python modules whose versions should be recorded in the simulation metadata. By default, all of the modules recorded by pypfilt.summary.metadata are included, as is the epifx package itself. **kwargs – Extra arguments to pass to pypfilt.summary.HDF5.

class epifx.summary.PrOutbreak(name='pr_epi')¶

Record the daily outbreak probability, defined as the sum of the weights of all particles in which an outbreak has been seeded.

Parameters:	name – the name of the table in the output file.

class epifx.summary.PeakMonitor¶

Record epidemic peak forecasts, for use with other statistics.

peak_size = None¶

A dictionary that maps observation systems to the size of each particle’s peak with respect to that system: peak_size[(unit, period)].

Note that this is only valid for tables to inspect in the finished() method, and not in the add_rows() method.

peak_date = None¶

A dictionary that maps observation systems to the date of each particle’s peak with respect to that system: peak_date[(unit, period)].

Note that this is only valid for tables to inspect in the finished() method, and not in the add_rows() method.

peak_time = None¶

A dictionary that maps observation systems to the time of each particle’s peak with respect to that system, measured in (fractional) days from the start of the forecasting period: peak_time[(unit, period)].

Note that this is only valid for tables to inspect in the finished() method, and not in the add_rows() method.

peak_weight = None¶

A dictionary that maps observation systems to the weight of each particle at the time that its peak occurs: peak_weight[(unit, period)].

Note that this is only valid for tables to inspect in the finished() method, and not in the add_rows() method.

expected_obs = None¶

The expected observation for each particle for the duration of the current simulation window.

Note that this is only valid for tables to inspect in each call to add_rows(), and not in a call to finished().

days_to(date)¶

Convert a date to the (fractional) number of days from the start of the forecasting period.

class epifx.summary.PeakForecastEnsembles(peak_monitor, name='peak_ensemble')¶

Record the weighted ensemble of peak size and time predictions for each forecasting simulation.

Parameters:	peak_monitor – an instance of PeakMonitor. name – the name of the table in the output file.

class epifx.summary.PeakForecastCIs(peak_monitor, probs=None, name='peak_cints')¶

Record fixed-probability central credible intervals for the peak size and time predictions.

Parameters:	peak_monitor – an instance of PeakMonitor. probs – an array of probabilities that define the size of each central credible interval. The default value is numpy.uint8([0, 50, 90, 95, 99, 100]). name – the name of the table in the output file.

class epifx.summary.PeakSizeAccuracy(peak_monitor, name='peak_size_acc', toln=None)¶

Record the accuracy of the peak size predictions against multiple accuracy tolerances.

Parameters:	peak_monitor – an instance of PeakMonitor. name – the name of the table in the output file. toln – The accuracy thresholds for peak size predictions, expressed as percentages of the true size. The default is np.array([10, 20, 25, 33]).

class epifx.summary.PeakTimeAccuracy(peak_monitor, name='peak_time_acc', toln=None)¶

Record the accuracy of the peak time predictions against multiple accuracy tolerances.

Parameters:	peak_monitor – an instance of PeakMonitor. name – the name of the table in the output file. toln – The accuracy thresholds for peak time predictions, expressed as numbers of days. The default is np.array([7, 10, 14]).

class epifx.summary.ExpectedObs(peak_monitor, probs=None, name='forecasts')¶

Record fixed-probability central credible intervals for the expected observations.

Parameters:	peak_monitor – an instance of PeakMonitor. probs – an array of probabilities that define the size of each central credible interval. The default value is numpy.uint8([0, 50, 90, 95, 99, 100]). name – the name of the table in the output file.

class epifx.summary.ObsLikelihood(name='obs_llhd')¶

Record the likelihood of each observation according to each particle. Note that this table registers its record_obs_llhd method in the parameter dictionary so that it can obtain the observation likelihoods.

Parameters:	name – the name of the table in the output file.

Generating forecasts¶

import datetime
import epifx

# Simulation parameters
p_exp = 1.0 / 36
px_scale = 100
params = epifx.default_params(p_exp, px_scale)
epifx.init_prng(params, seed=42)

# Simulate from the 1st of May to the 31st of October, 2015.
start = datetime.datetime(2015, 5, 1)
until = start + datetime.timedelta(days=183)

# Load the observations and define the observation model parameters.
obs_list = ...
params['obs'][obs_unit] = {...}

# Generate weekly forecasts for the first 9 weeks (and the start date).
fs_dates = [start + datetime.timedelta(days=week * 7)
            for week in range(10)]

# Summary statistics and output file.
summary = epifx.summary.make(params, obs_list))
out = "forecasts.hdf5"

# Generate the forecasts and save them to disk.
pypfilt.forecast(params, start, until, [obs_list], fs_dates, summary, out)

Comparing forecast outputs¶

Output files can be compared for equality, which is useful for ensuring that different systems produce identical results.

epifx.cmp.files(path1, path2, verbose=True, examine=None)¶

Compare two HDF5 files for identical simulation outputs.

Parameters:	path1 – The filename of the first HDF5 file. path2 – The filename of the second HDF5 file. verbose – Whether to report successful matches. examine – The data groups to examine for equality; the default is to examine the simulation outputs ('/data') and ignore the associated metadata ('/meta').
Returns:	True if the files contain identical simulation outputs, otherwise False.

This functionality is also provided as a command-line script:

cmp-output --help
cmp-output file1.hdf5 file2.hdf5

Note that simulation outputs have been observed to differ depending on the installed versions of the NumPy and SciPy packages, due to changes in the random number generator.

Observations from the entire population¶

The PopnCounts class provides a generic observation model for relating disease incidence to count data where the denominator is assumed to be the population size \(N\) (i.e., the denominator is assumed constant and is either not known or is known to be the population size \(N\)).

This observation model assumes that the relationship between observed case counts \(y_t\) and disease incidence in particle \(x_t\) follow a negative binomial distribution with mean \(\mathbb{E}[y_t]\) and dispersion parameter \(k\).

\[\begin{split}\mathcal{L}(y_t \mid x_t) &\sim NB(\mathbb{E}[y_t], k)\end{split}\]\[\begin{split}\mathbb{E}[y_t] &= (1 - p_\mathrm{inf}) \cdot bg_\mathrm{obs} + p_\mathrm{inf} \cdot p_\mathrm{obs} \cdot N\end{split}\]\[\begin{split}\operatorname{Var}[y_t] &= \mathbb{E}[y_t] + \frac{\left(\mathbb{E}[y_t]\right)^2}{k}\end{split}\]

The observation model parameters comprise:

• The background observation rate \(bg_\mathrm{obs}\);
• The probability of observing an infected individual \(p_\mathrm{obs}\); and
• The dispersion parameter \(k\), which controls the relationship between the mean \((\mathbb{E}[y_t])\) and the variance; as \(k \to \infty\) the distribution approaches the Poisson, as \(k \to 0\) the distribution becomes more and more over-dispersed with respect to the Poisson.

Probability mass functions for the PopnCounts observation model with different values of the dispersion parameter \(k\), for the expected value \(\mathbb{E}[y_t] = 200\) (vertical dashed line).

class epifx.obs.PopnCounts(obs_unit, obs_period, pr_obs_ix=None)¶

Generic observation model for relating disease incidence to count data where the denominator is assumed or known to be the population size.

Parameters:	obs_unit – A descriptive name for the data. obs_period – The observation period (in days). pr_obs_ix – The index of the model parameter that defines the observation probability (\(p_\mathrm{obs}\)). By default, the value in the parameters dictionary is used.

expect(params, op, period, pr_inf, prev, curr)¶

Calculate the expected observation value \(\mathbb{E}[y_t]\) for every particle \(x_t\).

log_llhd(params, op, obs, pr_indiv, curr, hist)¶

Calculate the log-likelihood \(\mathcal{l}(y_t \mid x_t)\) for the observation \(y_t\) (obs) and every particle \(x_t\).

If it is known (or suspected) that the observed value will increase in the future — when obs['incomplete'] == True — then the log-likehood \(\mathcal{l}(y > y_t \mid x_t)\) is calculated instead (i.e., the log of the survival function).

If an upper bound to this increase is also known (or estimated) — when obs['upper_bound'] is defined — then the log-likelihood \(\mathcal{l}(y_u \ge y > y_t \mid x_t)\) is calculated instead.

define_params(params, bg_obs, pr_obs, disp)¶

Add observation model parameters to the simulation parameters.

Parameters:	bg_obs – The background signal in the data (\(bg_\mathrm{obs}\)). pr_obs – The probability of observing an infected individual (\(p_\mathrm{obs}\)). disp – The dispersion parameter (\(k\)).

from_file(filename, year=None, date_col=u'to', value_col=u'count')¶

Load count data from a space-delimited text file with column headers defined in the first line.

Parameters:	filename – The file to read. year – Only returns observations for a specific year. The default behaviour is to return all recorded observations. date_col – The name of the observation date column. value_col – The name of the observation value column.
Returns:	A list of observations, ordered as per the original file.

Observations from population samples¶

The SampleCounts class provides a generic observation model for relating disease incidence to count data where the denominator is reported and may vary (e.g., weekly counts of all patients and the number that presented with influenza-like illness), and where the background signal is not a fixed value but rather a fixed proportion.

For this observation model, each observed value \((y_t)\) is the fraction of all patients \((N_t)\) that were identified as cases \((c_t)\). This observation model assumes that the relationship between the observed fraction \(y_t\), the denominator \(N_t\), and disease incidence in particle \(x_t\) follows a Beta-binomial distribution with probability of success \(\mathbb{E}[y_t]\) and dispersion parameter \(k\).

\[\begin{split}\mathcal{L}(c_t \mid N_t, x_t) &\sim \operatorname{BetaBin}(N_t, \mathbb{E}[y_t], k)\end{split}\]\[\begin{split}\mathbb{E}[c_t] &= N_t \cdot \mathbb{E}[y_t]\end{split}\]\[\begin{split}\operatorname{Var}[c_t] &= N_t \cdot p \cdot (1 - p) \cdot \left[ 1 + \frac{N_t - 1}{k + 1} \right]\end{split}\]\[\begin{split}\mathbb{E}[y_t] &= bg_\mathrm{obs} + p_\mathrm{inf} \cdot \kappa_\mathrm{obs}\end{split}\]\[\begin{split}y_t &= \frac{c_t}{N_t}\end{split}\]

The shape parameters \(\alpha\) and \(\beta\) satisfy:

\[\begin{split}p &= \mathbb{E}[y_t] = \frac{\alpha}{\alpha + \beta}\end{split}\]\[\begin{split}k &= \alpha + \beta\end{split}\]\[\begin{split}\implies \alpha &= p \cdot k\end{split}\]\[\begin{split}\implies \beta &= (1 - p) \cdot k\end{split}\]

The observation model parameters comprise:

• The background case fraction \(bg_\mathrm{obs}\);
• The slope \(\kappa_\mathrm{obs}\) of the relationship between disease incidence and the proportion of cases; and
• The dispersion parameter \(k\), which controls the relationship between the mean \((\mathbb{E}[y_t])\) and the variance; as \(k \to \infty\) the distribution approaches the binomial, as \(k \to 0\) the distribution becomes more and more over-dispersed with respect to the binomial.

Probability mass functions for the SampleCounts observation model with different values of the dispersion parameter \(k\), for the expected value \(\mathbb{E}[c_t] = 200\) (vertical dashed line) where \(\mathbb{E}[y_t] = 0.02\) and \(N_t = 10,\!000\).

class epifx.obs.SampleCounts(obs_unit, obs_period, k_obs_ix=None)¶

Generic observation model for relating disease incidence to count data where the denominator is reported.

Parameters:	obs_unit – A descriptive name for the data. obs_period – The observation period (in days). k_obs_ix – The index of the model parameter that defines the disease-related increase in observation rate (\(\kappa_\mathrm{obs}\)). By default, the value in the parameters dictionary is used.

expect(params, op, period, pr_inf, prev, curr)¶

Calculate the expected observation value \(\mathbb{E}[y_t]\) for every particle \(x_t\).

log_llhd(params, op, obs, pr_indiv, curr, hist)¶

Calculate the log-likelihood \(\mathcal{l}(y_t \mid x_t)\) for the observation \(y_t\) (obs) and every particle \(x_t\).

If it is known (or suspected) that the observed value will increase in the future — when obs['incomplete'] == True — then the log-likehood \(\mathcal{l}(y > y_t \mid x_t)\) is calculated instead (i.e., the log of the survival function).

If an upper bound to this increase is also known (or estimated) — when obs['upper_bound'] is defined — then the log-likelihood \(\mathcal{l}(y_u \ge y > y_t \mid x_t)\) is calculated instead.

define_params(params, bg_obs, k_obs, disp)¶

Add observation model parameters to the simulation parameters.

Parameters:	bg_obs – The background signal in the data (\(bg_\mathrm{obs}\)). k_obs – The increase in observation rate due to infected individuals (\(\kappa_\mathrm{obs}\)). disp – The dispersion parameter (\(k\)).

from_file(filename, year=None, date_col=u'to', value_col=u'cases', denom_col=u'patients')¶

Load count data from a space-delimited text file with column headers defined in the first line.

Note that returned observation values represent the fraction of patients that were counted as cases, not the absolute number of cases. The number of cases and the number of patients are recorded under the 'numerator' and 'denominator' keys, respectively.

Parameters:	filename – The file to read. year – Only returns observations for a specific year. The default behaviour is to return all recorded observations. date_col – The name of the observation date column. value_col – The name of the observation value column (reported as absolute values, not fractions). denom_col – The name of the observation denominator column.
Returns:	A list of observations, ordered as per the original file.
Raises:	ValueError – If a denominator or value is negative, or if the value exceeds the denominator.

Reading observations from disk¶

The from_file() methods provided by the observation models listed above use read_table() to read specific columns from data files. This function is intended to be a flexible wrapper around numpy.loadtxt, and should be sufficient for reading almost any whitespace-delimited data file with named columns.

epifx.obs.read_table(filename, columns, date_fmt=None, comment=u'#')¶

Load data from a space-delimited ASCII text file with column headers defined in the first non-comment line.

Parameters:	filename – The file to read. columns – The columns to read, represented as a list of (name, type) tuples; type should either be a built-in NumPy scalar type, or datetime.date or datetime.datetime (in which case string values will be automatically converted to datetime.datetime objects by datetime_conv()). date_fmt – A format string for parsing date columns; see datetime_conv() for details and the default format string. comment – Prefix string(s) that identify comment lines; either a single Unicode string or a tuple of Unicode strings.
Returns:	A NumPy structured array.
Raises:	ValueError – If the data file contains non-ASCII text (i.e., bytes greater than 127), because numpy.loadtxt cannot process non-ASCII data (e.g., see NumPy issues #3184, #4543, #4600, #4939).
Example:

The code below demonstrates how to read observations from a file that includes columns for the year, the observation date, the observed value, and free-text metadata (up to 20 characters in length).

import datetime
import numpy as np
import epifx.obs
columns = [('year', np.int32), ('date', datetime.datetime),
           ('count', np.int32), ('metadata', 'S20')]
df = epifx.obs.read_table('my-data-file.ssv', columns,
                          date_fmt='%Y-%m-%d')

epifx.obs.datetime_conv(text, fmt=u'%Y-%m-%d', encoding=u'utf-8')¶

Convert date strings to datetime.datetime instances. This is a convenience function intended for use with, e.g., numpy.loadtxt.

Parameters:	text – A string (bytes or Unicode) containing a date or date-time. fmt – A format string that defines the textual representation. See the Python strptime documentation for details. encoding – The byte string encoding (default is UTF-8).

0.4.3 (2016-09-16)¶

• This release requires pypfilt >= 0.4.3.
• Breaking change: the epifx.obs.SampleCounts observation model now uses a Beta-binomial distribution rather than a Beta distribution. Parameter names and definitions have been changed accordingly.
• Enhancement: consistently separate Unicode strings from bytes, and automatically convert NumPy field names into native strings.
• Enhancement: add support for incomplete data for which there may or may not be an upper bound (whether known in fact or estimated).
• Enhancement: record the likelihood of each observation according to each particle (see the epifx.summary.ObsLikelihood class).

0.4.2 (2016-06-16)¶

• Breaking change: replace the observation models added in epifx 0.4.1 with observations models for:
  • Count data where the denominator is known or assumed to be the entire population (epifx.obs.PopnCounts); and
  • Count data where the denominator is reported and may vary, and where the background signal is a fixed proportion (epifx.obs.SampleCounts).

0.4.1 (2016-04-22)¶

• Enhancement: provide generic negative binomial observation models for count data and for fractional/percentage data in epifx.obs.

0.4.0 (2016-04-22)¶

• This release requires pypfilt >= 0.4.0.

• Breaking change: models must define default parameter bounds by implementing the param_bounds method.

• Breaking change: model expectation functions now receive the previous and current state vectors, in addition to the infection probability vector. This means that expectation functions will need to change from:

  expect(params, unit, period, pr_inf)
  

  to:

  expect(params, unit, period, pr_inf, prev, curr)
  

• Enhancement: epifx.summary.make now passes additional keyword arguments to the pypfilt.summary.HDF5 constructor, allowing users to override default settings, such as first_day=False.

• Bug fix: ensure that infection probabilities are strictly non-negative.

• Bug fix: ensure that population invariants are enforced correctly.

• Bug fix: correctly scale the daily seeding probability.

• Add instructions for installing epifx in a virtual environment.

0.3.1 (2016-02-25)¶

• Bug fix: prevent a runtime error with params['epifx']['one_prng'] = True by correctly retrieving the pypfilt PRNG (params['resample']['rnd']).

0.3.0 (2016-02-23)¶

• This release requires pypfilt >= 0.3.0.
• Provide each summary statistic as a separate class.
• Inherit from the pypfilt simulation model base class.
• Host the documentation at Read The Docs.

0.2.0 (2015-11-17)¶

• Use an independent PRNG instance for model stochasticity, as distinct from the pypfilt PRNG instance (used for resampling). Note that this breaks backwards compatibility (in terms of producing identical outputs) and can be overridden by setting params['epifx']['one_prng'] = True.
• This release requires pypfilt >= 0.2.0.

0.1.10 (2015-07-09)¶

• Fix a bug where temporal forcing could cause a negative force of infection.

0.1.9 (2015-07-06)¶

• Add support for temporal forcing, modulated by the new parameter sigma.

0.1.8 (2015-06-18)¶

• Update the model parameter invariants for alpha, based on the priors for R0 and gamma.

0.1.7 (2015-06-18)¶

• Sample R0 and calculate alpha, rather than sampling alpha directly.

0.1.6 (2015-06-08)¶

• Avoid error messages if no logging handler is configured by the application.
• Default to comparing only the simulation outputs and ignore the metadata; this can be overridden by the --meta-data (-m) option.
• Build a universal wheel via python setup.py bdist_wheel, which supports both Python 2 and Python 3.
• This release requires pypfilt >= 0.1.2.

0.1.5 (2015-06-04)¶

• Record credible intervals for state variables (/data/state).

0.1.4 (2015-06-03)¶

• Reduce the minimum latent period to half a day.
• No longer require the simulation period to be included in the simulation parameters dictionary.
• Hide stderr output from spawned processes when obtaining git metadata, since the error messages have no relevance to the user (they only serve to indicate that the working directory is not part of a git repository).

0.1.3 (2015-06-01)¶

• Obtain git metadata from the working directory, if it is contained within a repository. This requires pypfilt >= 0.1.1.

  Note that sensible metadata will only be obtained if the working directory is not manipulated by program code (e.g., os.chdir).

0.1.2 (2015-06-01)¶

• Record the enforced limits on model parameters, so that output files include sufficient information to independently produce identical results.

  The default limits on beta and gamma are now identical to the domain of their default priors (1 to 3 days).

• Added options --verbose and --data-only to the cmp-output script.

• Ignore the command line (/meta/sim/cmdline) when comparing output files.

0.1.1 (2015-05-29)¶

• Added a script (cmp-output) that compares output files for identical simulation outputs and metadata.

0.1.0 (2015-05-29)¶

• Initial release.