Rt Estimation

Functions for estimating reproduction numbers. More...

Collaboration diagram for Rt Estimation:

Subgroups
	Estimates Smoothing
	Functions for smoothing estimates using Gaussian processes.
	Observation Model
	Functions for modeling the observation process.

Functions
vector	calculate_growth (vector infections, int seeding_time)
vector	calculate_Rt (vector infections, int seeding_time, vector gt_rev_pmf, int smooth)
real	R_to_r (real R, vector gt_rev_pmf, real abs_tol)
real	R_to_r_newton_step (real R, real r, vector pmf)
void	rt_lp (array[] real initial_infections_scale, vector bp_effects, array[] real bp_sd, int bp_n, array[] int cases, real initial_infections_guess)
vector	update_Rt (int t, real R0, vector noise, array[] int bps, vector bp_effects, int stationary)
	Update a vector of effective reproduction numbers (Rt) based on an intercept, breakpoints (i.e. a random walk), and a Gaussian process.

Description

Functions for estimating reproduction numbers.

This group contains functions for estimating and processing reproduction numbers (Rt), including updating Rt values, calculating growth rates, and converting between reproduction numbers and growth rates. Includes functions from rt.stan and related generated quantities.

Function Documentation

calculate_growth()

vector calculate_growth	(	vector	infections,
		int	seeding_time )

#include </github/workspace/inst/stan/functions/generated_quantities.stan>

Calculate growth rate

This function calculates the growth rate from a time series of infections by taking the log difference between consecutive time points.

Parameters

infections	Vector of infection counts
seeding_time	Number of time steps used for seeding

Returns

A vector of growth rates

Definition at line 63 of file generated_quantities.stan.

                                                             {
  int t = num_elements(infections);
  int ot = t - seeding_time;
  vector[t] log_inf = log(infections);
  vector[ot] growth =
    log_inf[(seeding_time + 1):t] - log_inf[seeding_time:(t - 1)];
  return(growth);
}

calculate_Rt()

vector calculate_Rt	(	vector	infections,
		int	seeding_time,
		vector	gt_rev_pmf,
		int	smooth )

#include </github/workspace/inst/stan/functions/generated_quantities.stan>

Generated Quantities Functions

Functions for calculating additional quantities from model outputs. Calculate Rt directly from inferred infections

This function estimates the reproduction number (Rt) using the Cori et al. approach, directly from a time series of infections. Optionally applies smoothing.

Parameters

infections	Vector of infection counts
seeding_time	Number of time steps used for seeding
gt_rev_pmf	Vector of reversed generation time PMF
smooth	Number of time steps to use for smoothing (0 for no smoothing)

Returns

A vector of reproduction numbers (Rt)

Definition at line 21 of file generated_quantities.stan.

                                                   {
  int t = num_elements(infections);
  int ot = t - seeding_time;
  vector[ot] R;
  vector[ot] sR;
  vector[ot] infectiousness = rep_vector(1e-5, ot);
  // calculate Rt using Cori et al. approach
  for (s in 1:ot) {
    infectiousness[s] += update_infectiousness(
      infections, gt_rev_pmf, seeding_time, s
    );
    R[s] = infections[s + seeding_time] / infectiousness[s];
  }
  if (smooth) {
    for (s in 1:ot) {
      real window = 0;
      sR[s] = 0;
      for (i in max(1, s - smooth):min(ot, s + smooth)) {
        sR[s] += R[i];
        window += 1;
      }
      sR[s] = sR[s] / window;
    }
  } else{
    sR = R;
  }
  return(sR);
}

References update_infectiousness().

Here is the call graph for this function:

R_to_r()

real R_to_r	(	real	R,
		vector	gt_rev_pmf,
		real	abs_tol )

#include </github/workspace/inst/stan/functions/rt.stan>

Estimate the growth rate r from reproduction number R

This function uses the Newton method to solve for the growth rate r that corresponds to a given reproduction number R, using the generation time distribution.

Code is based on Julia code from https://github.com/CDCgov/Rt-without-renewal/blob/d6344cc6e451e3e6c4188e4984247f890ae60795/EpiAware/test/predictive_checking/fast_approx_for_r.jl under Apache license 2.0.

Parameters

R	Reproduction number
gt_rev_pmf	Reversed probability mass function of the generation time
abs_tol	Absolute tolerance for the Newton solver

Returns

The estimated growth rate r

Definition at line 128 of file rt.stan.

                                                     {
  int gt_len = num_elements(gt_rev_pmf);
  vector[gt_len] gt_pmf = reverse(gt_rev_pmf);
  real mean_gt = dot_product(gt_pmf, linspaced_vector(gt_len, 0, gt_len - 1));
  real r = fmax((R - 1) / (R * mean_gt), -1);
  real step = abs_tol + 1;
  while (abs(step) > abs_tol) {
    step = R_to_r_newton_step(R, r, gt_pmf);
    r -= step;
  }
 
  return(r);
}

References R_to_r_newton_step().

Referenced by generate_infections().

Here is the call graph for this function:

Here is the caller graph for this function:

R_to_r_newton_step()

real R_to_r_newton_step	(	real	R,
		real	r,
		vector	pmf )

#include </github/workspace/inst/stan/functions/rt.stan>

Helper function for calculating r from R using Newton's method

This function performs a single Newton step in the iterative calculation of the growth rate r from the reproduction number R.

Code is based on Julia code from https://github.com/CDCgov/Rt-without-renewal/blob/d6344cc6e451e3e6c4188e4984247f890ae60795/EpiAware/test/predictive_checking/fast_approx_for_r.jl under Apache license 2.0.

Parameters

R	Reproduction number
r	Current estimate of the growth rate
pmf	Generation time probability mass function (first index: 0)

Returns

The Newton step for updating r

Definition at line 101 of file rt.stan.

                                                    {
  int len = num_elements(pmf);
  vector[len] zero_series = linspaced_vector(len, 0, len - 1);
  vector[len] exp_r = exp(-r * zero_series);
  real ret = (R * dot_product(pmf, exp_r) - 1) /
    (- R * dot_product(pmf .* zero_series, exp_r));
  return(ret);
}

Referenced by R_to_r().

Here is the caller graph for this function:

rt_lp()

void rt_lp	(	array[]real	initial_infections_scale,
		vector	bp_effects,
		array[]real	bp_sd,
		int	bp_n,
		array[]int	cases,
		real	initial_infections_guess )

#include </github/workspace/inst/stan/functions/rt.stan>

Calculate the log-probability of the reproduction number (Rt) priors

This function adds the log density contributions from priors on initial infections and breakpoint effects to the target.

Parameters

initial_infections_scale	Array of initial infection values
bp_effects	Vector of breakpoint effects
bp_sd	Array of breakpoint standard deviations
bp_n	Number of breakpoints
cases	Array of observed case counts
initial_infections_guess	Initial guess for infections based on cases

Definition at line 73 of file rt.stan.

                                          {
  //breakpoint effects on Rt
  if (bp_n > 0) {
    bp_sd[1] ~ normal(0, 0.1) T[0,];
    bp_effects ~ normal(0, bp_sd[1]);
  }
  initial_infections_scale ~ normal(initial_infections_guess, 2);
}

update_Rt()

vector update_Rt	(	int	t,
		real	R0,
		vector	noise,
		array[]int	bps,
		vector	bp_effects,
		int	stationary )

#include </github/workspace/inst/stan/functions/rt.stan>

Update a vector of effective reproduction numbers (Rt) based on an intercept, breakpoints (i.e. a random walk), and a Gaussian process.

Reproduction Number (Rt) Functions

This group of functions handles the calculation, updating, and conversion of reproduction numbers in the model. The reproduction number represents the average number of secondary infections caused by a single infected individual.

Parameters

t	Length of the time series
R0	Initial reproduction number
noise	Vector of Gaussian process noise values
bps	Array of breakpoint indices
bp_effects	Vector of breakpoint effects
stationary	Whether the Gaussian process is stationary (1) or non-stationary (0)

Returns

A vector of length t containing the updated Rt values

Definition at line 25 of file rt.stan.

                                                    {
  // define control parameters
  int bp_n = num_elements(bp_effects);
  int gp_n = num_elements(noise);
  // initialise intercept
  vector[t] logR = rep_vector(log(R0), t);
  //initialise breakpoints + rw
  if (bp_n) {
    vector[bp_n + 1] bp0;
    bp0[1] = 0;
    bp0[2:(bp_n + 1)] = cumulative_sum(bp_effects);
    logR = logR + bp0[bps];
  }
  //initialise gaussian process
  if (gp_n) {
    vector[t] gp = rep_vector(0, t);
    if (stationary) {
      gp[1:gp_n] = noise;
      // fix future gp based on last estimated
      if (t > gp_n) {
        gp[(gp_n + 1):t] = rep_vector(noise[gp_n], t - gp_n);
      }
    } else {
      gp[2:(gp_n + 1)] = noise;
      gp = cumulative_sum(gp);
    }
    logR = logR + gp;
  }
 
  return exp(logR);
}