In estimate_secondary(), does the obs refer to the % observed secondary cases among primary cases?

[avallecam]

About the obs argument:

• When used in estimate_infections(), the obs refers to the % observed reported cases among infections (i.e., under-report as used in tutorial episode on forecasting).
• When used in estimate_secondary(),
  • (A) Does the obs refer to the % observed secondary cases among primary cases? or
  • (B) Does the obs always refer to the % of observed primary observations (among infections)?

My mental map currently understands that:

Following this reference figure from Gosnic et al.:

infections --> %obs1 --> primary reports --> %obs2 --> secondary reports

• We can scale estimated infections by defining a known %obs1, given the observed primary cases.
• If we assume that secondary reports are deaths, then the %obs2 is analogous to the CFR
• As in the reprex below, we can get an estimate of the %obs2 by defining a scale in the observation model (obs) to get an estimate close to the one defined when secondary reports were generated by convolve_and_scale() of scaling := 0.4.

My questions are:

• Is the option (A) above an appropriate interpretation?

Reprex using Example in function documentation:

options(mc.cores = ifelse(interactive(), 4, 1))

# load data.table for manipulation
library(EpiNow2)
#> 
#> Attaching package: 'EpiNow2'
#> The following object is masked from 'package:stats':
#> 
#>     Gamma
library(data.table)
library(tidyverse)

#### Incidence data example ####

# make some example secondary incidence data
cases <- example_confirmed
cases <- as.data.table(cases)[, primary := confirm]
# Assume that only 40 percent of cases are reported
cases[, scaling := 0.4]
# Parameters of the assumed log normal delay distribution
cases[, meanlog := 1.8][, sdlog := 0.5]

# Simulate secondary cases
cases <- convolve_and_scale(cases, type = "incidence")
# print 
cases
#>            date confirm primary scaling meanlog sdlog index scaled       conv
#>          <Date>   <num>   <num>   <num>   <num> <num> <int>  <num>      <num>
#>   1: 2020-02-22      14      14     0.4     1.8   0.5     1    5.6   5.600000
#>   2: 2020-02-23      62      62     0.4     1.8   0.5     2   24.8   5.827560
#>   3: 2020-02-24      53      53     0.4     1.8   0.5     3   21.2   8.801043
#>   4: 2020-02-25      97      97     0.4     1.8   0.5     4   38.8  12.938343
#>   5: 2020-02-26      93      93     0.4     1.8   0.5     5   37.2  16.590082
#>  ---                                                                         
#> 126: 2020-06-26     296     296     0.4     1.8   0.5   126  118.4  91.432820
#> 127: 2020-06-27     255     255     0.4     1.8   0.5   127  102.0  95.786381
#> 128: 2020-06-28     175     175     0.4     1.8   0.5   128   70.0 103.811968
#> 129: 2020-06-29     174     174     0.4     1.8   0.5   129   69.6 109.093575
#> 130: 2020-06-30     126     126     0.4     1.8   0.5   130   50.4 109.029626
#>      secondary
#>          <int>
#>   1:         5
#>   2:         5
#>   3:         8
#>   4:        12
#>   5:        16
#>  ---          
#> 126:        91
#> 127:        95
#> 128:       103
#> 129:       109
#> 130:       109


# Fig. A: With DEFAULT scale in observation model Fixed(1) -----------------------------------

# fit model to example data specifying a weak prior for fraction reported
# with a secondary case
inc <- estimate_secondary(
  cases[1:60],
  obs = obs_opts(
    # scale = Normal(mean = 0.2, sd = 0.2), # Fixed(1) default 
    week_effect = FALSE
  )
)
#> WARN [2025-06-18 16:58:43] estimate_secondary (chain: 1): There were 2 divergent transitions after warmup. See
#> https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
#> to find out why this is a problem and how to eliminate them. - 
#> WARN [2025-06-18 16:58:44] estimate_secondary (chain: 1): Examine the pairs() plot to diagnose sampling problems
#>  -
# plot(inc, primary = TRUE)
frac_obs_var <- paste0("params[", inc$data$frac_obs_id, "]")
inc$posterior[variable == frac_obs_var]
#>     variable      mean     se_mean         sd lower_90  lower_50  lower_20
#>       <char>     <num>       <num>      <num>    <num>     <num>     <num>
#> 1: params[1] 0.5506713 0.001577778 0.05142625 0.472063 0.5155377 0.5340829
#>       median  upper_20  upper_50  upper_90
#>        <num>     <num>     <num>     <num>
#> 1: 0.5470943 0.5620268 0.5827512 0.6412143

# Fig. B: With DEFINED scale in observation model Normal(0.2, 0.2) --------------------------------------

# fit model to example data specifying a weak prior for fraction reported
# with a secondary case
inc <- estimate_secondary(
  cases[1:60],
  obs = obs_opts(
    scale = Normal(mean = 0.2, sd = 0.2), # Defined scale
    week_effect = FALSE
  )
)
# plot(inc, primary = TRUE)
frac_obs_var <- paste0("params[", inc$data$frac_obs_id, "]")
inc$posterior[variable == frac_obs_var]
#>     variable      mean      se_mean          sd  lower_90  lower_50  lower_20
#>       <char>     <num>        <num>       <num>     <num>     <num>     <num>
#> 1: params[1] 0.4001165 4.475948e-05 0.001851689 0.3970677 0.3988574 0.3996644
#>       median  upper_20  upper_50  upper_90
#>        <num>     <num>     <num>     <num>
#> 1: 0.4001262 0.4005981 0.4013479 0.4032705

^{Created on 2025-06-18 with reprex v2.1.1}

Fig. A: With default scale in observation model (Fixed(1))

plot(inc, primary = TRUE)

Captures the delay but not the scale. (delay defined as default in estimate_secondary())

Fig. B: With defined scale in obs (Normal(mean = 0.2, sd = 0.2))

plot(inc, primary = TRUE)

Captures the delay and the scale.

[jamesmbaazam]

Hi @avallecam Sorry I missed this question.

The obs$scale argument accounts for the level of under-reporting in the primary observations before being used to estimate the secondary outcomes.