EpiNow2 Stan Functions
generated_quantities.stan
Go to the documentation of this file.
1/**
2 * Generated Quantities Functions
3 *
4 * Functions for calculating additional quantities from model outputs.
5 */
6
7/**
8 * Calculate Rt directly from inferred infections
9 *
10 * This function estimates the reproduction number (Rt) using the Cori et al. approach,
11 * directly from a time series of infections. Optionally applies smoothing.
12 *
13 * @param infections Vector of infection counts
14 * @param seeding_time Number of time steps used for seeding
15 * @param gt_rev_pmf Vector of reversed generation time PMF
16 * @param smooth Number of time steps to use for smoothing (0 for no smoothing)
17 * @return A vector of reproduction numbers (Rt)
18 *
19 * @ingroup rt_estimation
20 */
21vector calculate_Rt(vector infections, int seeding_time,
22 vector gt_rev_pmf, int smooth) {
23 int t = num_elements(infections);
24 int ot = t - seeding_time;
25 vector[ot] R;
26 vector[ot] sR;
27 vector[ot] infectiousness = rep_vector(1e-5, ot);
28 // calculate Rt using Cori et al. approach
29 for (s in 1:ot) {
30 infectiousness[s] += update_infectiousness(
31 infections, gt_rev_pmf, seeding_time, s
32 );
33 R[s] = infections[s + seeding_time] / infectiousness[s];
34 }
35 if (smooth) {
36 for (s in 1:ot) {
37 real window = 0;
38 sR[s] = 0;
39 for (i in max(1, s - smooth):min(ot, s + smooth)) {
40 sR[s] += R[i];
41 window += 1;
42 }
43 sR[s] = sR[s] / window;
44 }
45 } else{
46 sR = R;
47 }
48 return(sR);
49}
50
51/**
52 * Calculate growth rate
53 *
54 * This function calculates the growth rate from a time series of infections
55 * by taking the log difference between consecutive time points.
56 *
57 * @param infections Vector of infection counts
58 * @param seeding_time Number of time steps used for seeding
59 * @param gt_rev_pmf Vector of reversed generation time PMF
60 * @param growth_method Either 0 (log derivative of new infections) or 1 (log
61 * derivative of infectiousness, see Parag et al. 2022)
62 *
63 * @return A vector of growth rates
64 *
65 * @ingroup rt_estimation
66 */
67vector calculate_growth(vector infections, int seeding_time,
68 vector gt_rev_pmf, int growth_method) {
69 if (growth_method == 0) {
70 return(calculate_growth_infections(infections, seeding_time));
71 } else if (growth_method == 1) {
72 return(calculate_growth_infness(infections, seeding_time, gt_rev_pmf));
73 } else {
74 reject("growth_method must be 0 (infections) or 1 (infectiousness).");
75 }
76}
77
78/**
79 * Calculate growth rate
80 *
81 * This function calculates the growth rate from a time series of infections
82 * by taking the log difference between consecutive time points.
83 *
84 * @param infections Vector of infection counts
85 * @param seeding_time Number of time steps used for seeding
86 *
87 * @return A vector of growth rates
88 *
89 * @ingroup rt_estimation
90 */
91vector calculate_growth_infections(vector infections, int seeding_time) {
92 int t = num_elements(infections);
93 int ot = t - seeding_time;
94 int start = 1 + seeding_time;
95 vector[t] log_inf = log(infections);
96 vector[ot - 1] growth = log_inf[(1+start):t] - log_inf[start:(t - 1)];
97 return(growth);
98}
99
100/**
101 * Calculate growth rate using approach by Parag et al. 2022
102 * (https://doi.org/10.1111/rssa.12867)
103 *
104 * This function calculates the growth rate from a time series of infections
105 * by taking the log derivative on the infectiousness and shifting it by the
106 * mean generation time.
107 *
108 * @param infections Vector of infection counts
109 * @param seeding_time Number of time steps used for seeding
110 * @param gt_rev_pmf Vector of reversed generation time PMF
111 *
112 * @return A vector of growth rates
113 *
114 * @ingroup rt_estimation
115 */
116vector calculate_growth_infness(vector infections, int seeding_time,
117 vector gt_rev_pmf) {
118 int t = num_elements(infections);
119 int ot = t - seeding_time;
120 int start = 1 + seeding_time;
121 if (ot <= 1) {
122 reject("seeding_time must >1 time step shorter than infections vector.");
123 }
124 // infectiousness
125 vector[ot] infness_log = rep_vector(1e-5, ot);
126 for (s in 1:ot) {
127 infness_log[s] += log(update_infectiousness(
128 infections, gt_rev_pmf, seeding_time, s
129 ));
130 }
131 // mean generation time, will always be >= 1
132 int gt_length = num_elements(gt_rev_pmf);
133 int mean_gen = to_next_int_index( // round weighted mean to next int
134 dot_product(reverse(linspaced_vector(gt_length, 1, gt_length)), gt_rev_pmf)
135 );
136 // growth rate
137 vector[ot - 1] growth = infness_log[2:ot] - infness_log[1:(ot - 1)];
138 // shift by mean_gen (most recent growth rates remain undefined)
139 growth[1:(ot - 1 - mean_gen)] = growth[(1 + mean_gen):(ot - 1)];
140 growth[(ot - mean_gen):(ot - 1)] = rep_vector(not_a_number(), mean_gen);
141 return(growth);
142}
real update_infectiousness(vector infections, vector gt_rev_pmf, int seeding_time, int index)
vector calculate_Rt(vector infections, int seeding_time, vector gt_rev_pmf, int smooth)
vector calculate_growth_infections(vector infections, int seeding_time)
vector calculate_growth(vector infections, int seeding_time, vector gt_rev_pmf, int growth_method)
vector calculate_growth_infness(vector infections, int seeding_time, vector gt_rev_pmf)
int to_next_int_index(real x)
Definition helpers.stan:1