Observation Model

Functions for modeling the observation process. More...

Collaboration diagram for Observation Model:

Functions
vector	accumulate_reports (vector reports, array[] int accumulate)
vector	day_of_week_effect (vector reports, array[] int day_of_week, vector effect)
vector	report_log_lik (array[] int cases, vector reports, real dispersion, int model_type, real weight)
void	report_lp (array[] int cases, array[] int case_times, vector reports, real dispersion, int model_type, real weight)
array[] int	report_rng (vector reports, real dispersion, int model_type)
vector	scale_obs (vector reports, real frac_obs)
vector	truncate_obs (vector reports, vector trunc_rev_cmf, int reconstruct)
void	truncation_lp (array[] real truncation_mean, array[] real truncation_sd, array[] real trunc_mean_mean, array[] real trunc_mean_sd, array[] real trunc_sd_mean, array[] real trunc_sd_sd)

Description

Functions for modeling the observation process.

This group contains functions related to how infections or reports are observed, including day-of-week effects, reporting delays, and various likelihood models.

Function Documentation

accumulate_reports()

vector accumulate_reports	(	vector	reports,
		array[]int	accumulate )

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

Accumulate reports according to a binary flag at each time point

This function accumulates reports according to a binary flag at each time point.

Parameters

reports	Vector of expected reports.
accumulate	Array of integers indicating, for each time point, whether to accumulate or not.

Returns

A vector of accumulated reports.

Definition at line 161 of file observation_model.stan.

                                                                  {
  int ot_h = num_elements(reports); // number of reporting time points modelled
  vector[ot_h] accumulated_reports = reports;
  for (i in 1:(ot_h - 1)) {
    if (accumulate[i]) { // first observation gets ignored when accumulating
      accumulated_reports[i + 1] += accumulated_reports[i];
    }
  }
  return accumulated_reports;
}

day_of_week_effect()

vector day_of_week_effect	(	vector	reports,
		array[]int	day_of_week,
		vector	effect )

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

Apply day of the week effect to reports

This function applies a day of the week effect to a vector of reports.

Parameters

reports	Vector of reports to be adjusted.
day_of_week	Array of integers representing the day of the week for each report.
effect	Vector of day of week effects.

Returns

A vector of reports adjusted for day of the week effects.

Definition at line 15 of file observation_model.stan.

                                         {
  int wl = num_elements(effect);
  vector[wl] scaled_effect = wl * effect;
  return reports .* scaled_effect[day_of_week];
}

report_log_lik()

vector report_log_lik	(	array[]int	cases,
		vector	reports,
		real	dispersion,
		int	model_type,
		real	weight )

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

Calculate log likelihood for reported cases

This function calculates the log likelihood for reported cases based on the specified model type.

Parameters

cases	Array of integer observed cases.
reports	Vector of expected reports.
dispersion	Array of real values for reporting overdispersion.
model_type	Integer indicating the model type (0 for Poisson, >0 for Negative Binomial).
weight	Real value for weighting the log likelihood contribution.

Returns

A vector of log likelihoods for each time point.

Definition at line 189 of file observation_model.stan.

                                                                    {
  int t = num_elements(reports);
  vector[t] log_lik;
 
  // defer to poisson if phi is large, to avoid overflow
  if (model_type == 0) {
    for (i in 1:t) {
      log_lik[i] = poisson_lpmf(cases[i] | reports[i]) * weight;
    }
  } else {
    real phi = inv_square(dispersion);
    for (i in 1:t) {
      log_lik[i] = neg_binomial_2_lpmf(
        cases[i] | reports[i], dispersion
      ) * weight;
    }
  }
  return(log_lik);
}

report_lp()

void report_lp	(	array[]int	cases,
		array[]int	case_times,
		vector	reports,
		real	dispersion,
		int	model_type,
		real	weight )

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

Update log density for reported cases

This function updates the log density for reported cases based on the specified model type.

Parameters

cases	Array of integer observed cases.
case_times	Array of integer time indices for observed cases.
reports	Vector of expected reports.
dispersion	Real values for reporting overdispersion.
model_type	Integer indicating the model type (0 for Poisson, >0 for Negative Binomial).
weight	Real value for weighting the log density contribution.

Definition at line 125 of file observation_model.stan.

                                                             {
  int n = num_elements(case_times); // number of observations
  vector[n] obs_reports = reports[case_times]; // reports at observation time
  if (model_type) {
    real phi = inv_square(dispersion);
    if (weight == 1) {
      cases ~ neg_binomial_2(obs_reports, phi);
    } else {
      target += neg_binomial_2_lpmf(
        cases | obs_reports, phi
      ) * weight;
    }
  } else {
    if (weight == 1) {
      cases ~ poisson(obs_reports);
    } else {
      target += poisson_lpmf(cases | obs_reports) * weight;
    }
  }
}

report_rng()

array[] int report_rng	(	vector	reports,
		real	dispersion,
		int	model_type )

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

Generate random samples of reported cases

This function generates random samples of reported cases based on the specified model type.

Parameters

reports	Vector of expected reports.
dispersion	Real value for reporting overdispersion.
model_type	Integer indicating the model type (0 for Poisson, >0 for Negative Binomial).

Returns

An array of integer sampled reports.

Definition at line 257 of file observation_model.stan.

                                                                        {
  int t = num_elements(reports);
  array[t] int sampled_reports;
  real phi = 1e5;
  if (model_type) {
    phi = inv_square(dispersion);
  }
 
  for (s in 1:t) {
    sampled_reports[s] = neg_binomial_2_safe_rng(reports[s], phi);
  }
  return(sampled_reports);
}

References neg_binomial_2_safe_rng().

Here is the call graph for this function:

scale_obs()

vector scale_obs	(	vector	reports,
		real	frac_obs )

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

Scale observations by fraction reported

This function scales a vector of reports by a fraction observed.

Parameters

reports	Vector of reports to be scaled.
frac_obs	Real value representing the fraction observed.

Returns

A vector of scaled reports.

Definition at line 34 of file observation_model.stan.

                                                {
  int t = num_elements(reports);
  vector[t] scaled_reports;
  scaled_reports = reports * frac_obs;
  return(scaled_reports);
}

truncate_obs()

vector truncate_obs	(	vector	reports,
		vector	trunc_rev_cmf,
		int	reconstruct )

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

Truncate observed data by a truncation distribution

This function truncates a vector of reports based on a truncation distribution.

Parameters

reports	Vector of reports to be truncated.
trunc_rev_cmf	Vector representing the reverse cumulative mass function of the truncation distribution.
reconstruct	Integer flag indicating whether to reconstruct (1) or truncate (0) the data.

Returns

A vector of truncated reports.

Definition at line 57 of file observation_model.stan.

                                                                           {
  int t = num_elements(reports);
  int trunc_max = num_elements(trunc_rev_cmf);
  vector[t] trunc_reports = reports;
  // Calculate cmf of truncation delay
  int joint_max = min(t, trunc_max);
  int first_t = t - joint_max + 1;
  int first_trunc = trunc_max - joint_max + 1;
 
  // Apply cdf of truncation delay to truncation max last entries in reports
  if (reconstruct) {
    trunc_reports[first_t:t] ./= trunc_rev_cmf[first_trunc:trunc_max];
  } else {
    trunc_reports[first_t:t] .*= trunc_rev_cmf[first_trunc:trunc_max];
  }
  return(trunc_reports);
}

truncation_lp()

void truncation_lp	(	array[]real	truncation_mean,
		array[]real	truncation_sd,
		array[]real	trunc_mean_mean,
		array[]real	trunc_mean_sd,
		array[]real	trunc_sd_mean,
		array[]real	trunc_sd_sd )

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

Update log density for truncation distribution priors

This function updates the log density for truncation distribution priors.

Parameters

truncation_mean	Array of real values for truncation mean.
truncation_sd	Array of real values for truncation standard deviation.
trunc_mean_mean	Array of real values for mean of truncation mean prior.
trunc_mean_sd	Array of real values for standard deviation of truncation mean prior.
trunc_sd_mean	Array of real values for mean of truncation standard deviation prior.
trunc_sd_sd	Array of real values for standard deviation of truncation standard deviation prior.

Definition at line 93 of file observation_model.stan.

                                                                         {
  int truncation = num_elements(truncation_mean);
  if (truncation) {
    if (trunc_mean_sd[1] > 0) {
      // uncertain mean
      truncation_mean ~ normal(trunc_mean_mean, trunc_mean_sd);
    }
    if (trunc_sd_sd[1] > 0) {
      // uncertain sd
      truncation_sd ~ normal(trunc_sd_mean, trunc_sd_sd);
    }
  }
}