how to get the mean_sd and sd_sd from mean_ci_limits and sd_ci_limits for uncertain distributions in EpiNow2?

[avallecam]

For EpiNow2 we can use uncertain distributions by adding mean_sd and sd_sd as part of the dist_spec().

For mean_sd I followed the steps below matching cochran handbook using mean_ci_limits from an epidist object. Please, see the code below for your review.

How to get the sd_sd? I found a reference at stackoverflow. Can you help me get the code for this?

#' if
#' sem = standard error of the mean
#' sd = standard deviation
#' n = sample size
#' 1.96 = critical value for a significant level of 5% from qnorm(p = 0.975)
#' then
#' sem = sd / sqrt(n)
#' and
#' 95%ci = mean +- 1.96*sem
#' where
#' precision = 1.96*sem
#'
#' thus, if we have the 95% ci width (mean_ci_width)
#' mean_ci_width = 2 * precision
#' mean_ci_width = 2 * 1.96 * sem
#' mean_ci_width = 2 * 1.96 * sd / sqrt(n)
#' we have
#' sd = sqrt(n) * (mean_ci_width / (2 * 1.96))

covid_lnorm <- 
  epiparameter::epidist_db(
    disease = "covid",
    epi_dist = "serial",
    author = "Nishiura",
    single_epidist = T
  )
#> Using Nishiura H, Linton N, Akhmetzhanov A (2020). "Serial interval of novel
#> coronavirus (COVID-19) infections." _International Journal of
#> Infectious Diseases_. doi:10.1016/j.ijid.2020.02.060
#> <https://doi.org/10.1016/j.ijid.2020.02.060>.. 
#> To retrieve the short citation use the 'get_citation' function

covid_lnorm
#> Disease: COVID-19
#> Pathogen: SARS-CoV-2
#> Epi Distribution: serial interval
#> Study: Nishiura H, Linton N, Akhmetzhanov A (2020). "Serial interval of novel
#> coronavirus (COVID-19) infections." _International Journal of
#> Infectious Diseases_. doi:10.1016/j.ijid.2020.02.060
#> <https://doi.org/10.1016/j.ijid.2020.02.060>.
#> Distribution: lnorm
#> Parameters:
#>   meanlog: 1.386
#>   sdlog: 0.568

# mean
covid_lnorm$summary_stats$mean
#> [1] 4.7

# mean_ci width
covid_lnorm$summary_stats$mean_ci
#> [1] 95
covid_lnorm$summary_stats$mean_ci_limits
#> [1] 3.7 6.0
mean_ci_limits_num <- covid_lnorm$summary_stats$mean_ci_limits
mean_ci_width <- mean_ci_limits_num[2] - mean_ci_limits_num[1]
mean_ci_width
#> [1] 2.3

# from paper
covid_lnorm_sample <- covid_lnorm$metadata$sample_size 
stats::qt(p = 0.975,df = covid_lnorm_sample-1)
#> [1] 2.051831
stats::qt(p = 0.025,df = covid_lnorm_sample-1)
#> [1] -2.051831
t_095 <- stats::qt(p = 0.975,df = covid_lnorm_sample-1)

# mean_sd
covid_lnorm_mean_sd <- sqrt(covid_lnorm_sample) * (mean_ci_width / 2*t_095)
covid_lnorm_mean_sd
#> [1] 12.48586

^{Created on 2023-11-30 with reprex v2.0.2}

[joshwlambert]

@avallecam Thanks for question. I will need some time to think about the best approach.

Tagging @adamkucharski & @sbfnk for input as well.

[sbfnk]

I'm not sure there's a hard and fast method - in <EpiNow2> this represents a Bayesian prior which is assumed to be normally distributed with a given sd. If the CI represents e.g. frequentist confidence intervals then perhaps using the standard error directly would be more appropriate instead of dividing by sqrt(n)?

[avallecam]

Thank you for the clarification. I think I understand your point, we can look for the sem defined as the standard deviation of the sampling distribution of sample means.

Please, see below the re-written code:

# if
# sem = standard error of the mean = standard deviation of the sample means distribution
# sd = standard deviation of the sample distribution
# ci = confidence intervals of the population mean
# n = sample size
# 1.96 = critical value for a significant level of 5% from qnorm(p = 0.975)
# then
# 95%ci = mean +- 1.96*sem
# precision = 1.96*sem
# 
# for the mean_sd in {EpiNow2}
# we need the standard deviation of the mean distribution
# represented as a Bayesian prior which is assumed to be normally distributed with a given sd
# then,
# for the mean_sd we can use the sem
# thus,
# with mean_ci_width as the 95% ci width of the mean
# mean_ci_width = 2 * precision
# mean_ci_width = 2 * 1.96 * sem
# we have
# sem = (mean_ci_width / (2 * 1.96))

covid_lnorm <-
  epiparameter::epidist_db(
    disease = "covid",
    epi_dist = "serial",
    author = "Nishiura",
    single_epidist = T
  )
#> Using Nishiura H, Linton N, Akhmetzhanov A (2020). "Serial interval of novel
#> coronavirus (COVID-19) infections." _International Journal of
#> Infectious Diseases_. doi:10.1016/j.ijid.2020.02.060
#> <https://doi.org/10.1016/j.ijid.2020.02.060>.. 
#> To retrieve the short citation use the 'get_citation' function

covid_lnorm
#> Disease: COVID-19
#> Pathogen: SARS-CoV-2
#> Epi Distribution: serial interval
#> Study: Nishiura H, Linton N, Akhmetzhanov A (2020). "Serial interval of novel
#> coronavirus (COVID-19) infections." _International Journal of
#> Infectious Diseases_. doi:10.1016/j.ijid.2020.02.060
#> <https://doi.org/10.1016/j.ijid.2020.02.060>.
#> Distribution: lnorm
#> Parameters:
#>   meanlog: 1.386
#>   sdlog: 0.568

# mean
covid_lnorm$summary_stats$mean
#> [1] 4.7

# mean_ci width
covid_lnorm$summary_stats$mean_ci
#> [1] 95
covid_lnorm$summary_stats$mean_ci_limits
#> [1] 3.7 6.0
mean_ci_limits_num <- covid_lnorm$summary_stats$mean_ci_limits
mean_ci_width <- mean_ci_limits_num[2] - mean_ci_limits_num[1]
mean_ci_width
#> [1] 2.3

# from paper
covid_lnorm_sample <- covid_lnorm$metadata$sample_size
covid_lnorm_sample
#> [1] 28
# stats::qt(p = 0.975,df = covid_lnorm_sample-1)
# stats::qt(p = 0.025,df = covid_lnorm_sample-1)
t_095 <- stats::qt(p = 0.975,df = covid_lnorm_sample-1)

# mean_sd
covid_lnorm_mean_sd <- (mean_ci_width / 2*t_095)
covid_lnorm_mean_sd
#> [1] 2.359605

^{Created on 2023-12-16 with reprex v2.0.2}

This last mean_sd looks like a more fair value.

[avallecam]

Can we follow the same logic for the sd_sd?

# sd_sd
covid_lnorm$summary_stats$sd
#> [1] 2.9
sd_ci_limits_num <- covid_lnorm$summary_stats$sd_ci_limits
sd_ci_width <- sd_ci_limits_num[2] - sd_ci_limits_num[1]
sd_ci_width
#> [1] 3
covid_lnorm_sd_sd <- (sd_ci_width / 2*t_095)
covid_lnorm_sd_sd
#> [1] 3.077746

^{Created on 2023-12-16 with reprex v2.0.2}

[adamkucharski]

I also wonder how much difference including mean_sd>0 makes in practice, given the variance of the mean is much smaller than the variance of the delay distribution, and count data (like daily reported cases) sums over these delays? If we set mean_sd=0, rather than =0.1, intuitively I'd expect a very similar Rt curve (but may well be wrong!)

[avallecam]

@jamesmbaazam @joshwlambert in case related with a chat about uncertainty from epiparameter to epinow2

[avallecam]

Following epiforecasts/EpiNow2#504 and what worked in epiverse-trace/episoap#130 this may be how the current connection looks like

# howto epiparameter uncertain to epinow2 ---------------------------------

library(epiparameter)
library(EpiNow2)
#> 
#> Attaching package: 'EpiNow2'
#> The following object is masked from 'package:epiparameter':
#> 
#>     discretise
#> The following object is masked from 'package:stats':
#> 
#>     Gamma

# if
# sem = standard error of the mean = standard deviation of the sample means distribution
# sd = standard deviation of the sample distribution
# ci = confidence intervals of the population mean
# n = sample size
# 1.96 = critical value for a significant level of 5% from qnorm(p = 0.975)
# then
# 95%ci = mean +- 1.96*sem
# precision = 1.96*sem
# 
# for the mean_sd in {EpiNow2}
# we need the standard deviation of the mean distribution
# represented as a Bayesian prior which is assumed to be normally distributed with a given sd
# then,
# for the mean_sd we can use the sem
# thus,
# with mean_ci_width as the 95% ci width of the mean
# mean_ci_width = 2 * precision
# mean_ci_width = 2 * 1.96 * sem
# we have
# sem = (mean_ci_width / (2 * 1.96))

covid_lnorm <-
  epiparameter::epidist_db(
    disease = "covid",
    epi_dist = "serial",
    author = "Nishiura",
    single_epidist = T
  )
#> Using Nishiura H, Linton N, Akhmetzhanov A (2020). "Serial interval of novel
#> coronavirus (COVID-19) infections." _International Journal of
#> Infectious Diseases_. doi:10.1016/j.ijid.2020.02.060
#> <https://doi.org/10.1016/j.ijid.2020.02.060>.. 
#> To retrieve the short citation use the 'get_citation' function

covid_lnorm
#> Disease: COVID-19
#> Pathogen: SARS-CoV-2
#> Epi Distribution: serial interval
#> Study: Nishiura H, Linton N, Akhmetzhanov A (2020). "Serial interval of novel
#> coronavirus (COVID-19) infections." _International Journal of
#> Infectious Diseases_. doi:10.1016/j.ijid.2020.02.060
#> <https://doi.org/10.1016/j.ijid.2020.02.060>.
#> Distribution: lnorm
#> Parameters:
#>   meanlog: 1.386
#>   sdlog: 0.568

# mean
covid_lnorm$summary_stats$mean
#> [1] 4.7


# get the mean sd ---------------------------------------------------------

# mean_ci width
covid_lnorm$summary_stats$mean_ci
#> [1] 95
covid_lnorm$summary_stats$mean_ci_limits
#> [1] 3.7 6.0
mean_ci_limits_num <- covid_lnorm$summary_stats$mean_ci_limits
mean_ci_width <- mean_ci_limits_num[2] - mean_ci_limits_num[1]
mean_ci_width
#> [1] 2.3

# from paper
covid_lnorm_sample <- covid_lnorm$metadata$sample_size
covid_lnorm_sample
#> [1] 28
# stats::qt(p = 0.975,df = covid_lnorm_sample-1)
# stats::qt(p = 0.025,df = covid_lnorm_sample-1)
t_095 <- stats::qt(p = 0.975,df = covid_lnorm_sample-1)

# mean_sd
covid_lnorm_mean_sd <- (mean_ci_width / 2*t_095)
covid_lnorm_mean_sd
#> [1] 2.359605

# get the sd sd -----------------------------------------------------------

# sd_sd
covid_lnorm$summary_stats$sd
#> [1] 2.9
sd_ci_limits_num <- covid_lnorm$summary_stats$sd_ci_limits
sd_ci_width <- sd_ci_limits_num[2] - sd_ci_limits_num[1]
sd_ci_width
#> [1] 3
covid_lnorm_sd_sd <- (sd_ci_width / 2*t_095)
covid_lnorm_sd_sd
#> [1] 3.077746

# for epinow2 -------------------------------------------------------------

## discretise continuous distribution --------------------------------------

covid_lnorm_discrete <- 
  epiparameter::discretise(covid_lnorm)

covid_lnorm_discrete
#> Disease: COVID-19
#> Pathogen: SARS-CoV-2
#> Epi Distribution: serial interval
#> Study: Nishiura H, Linton N, Akhmetzhanov A (2020). "Serial interval of novel
#> coronavirus (COVID-19) infections." _International Journal of
#> Infectious Diseases_. doi:10.1016/j.ijid.2020.02.060
#> <https://doi.org/10.1016/j.ijid.2020.02.060>.
#> Distribution: discrete lnorm
#> Parameters:
#>   meanlog: 1.386
#>   sdlog: 0.568

## get maximum value from discrete distribution ----------------------------

covid_lnorm_discrete_max <- 
  covid_lnorm_discrete$prob_dist$q(p = 0.999)

covid_lnorm_discrete_max
#> [1] 23

## plug in to lognormal ----------------------------------------------------

LogNormal(
  mean = Normal(
    mean = covid_lnorm$summary_stats$mean,
    sd = covid_lnorm_mean_sd), 
  sd = Normal(
    mean = covid_lnorm$summary_stats$sd,
    sd = covid_lnorm_sd_sd),
  max = covid_lnorm_discrete_max
)
#> Warning in new_dist_spec(params, "lognormal"): Uncertain lognormal distribution
#> specified in terms of parameters that are not the "natural" parameters of the
#> distribution (meanlog, sdlog). Converting using a crude and very approximate
#> method that is likely to produce biased results. If possible, it is preferable
#> to specify the distribution directly in terms of the natural parameters.
#> - lognormal distribution (max: 23):
#>   meanlog:
#>     - normal distribution:
#>       mean:
#>         1.4
#>       sd:
#>         1.8
#>   sdlog:
#>     - normal distribution:
#>       mean:
#>         0.57
#>       sd:
#>         1.1

# a lognormal with max value is expected to run for epinow()

^{Created on 2024-02-13 with reprex v2.0.2}

[avallecam]

Is this making the same mistake as in epiverse-trace/episoap#130 (comment)?

[avallecam]

I think this is not, because this is using the summary statistics uncertainty (normally distributed) for the uncertainty of distribution parameters (lognormal distribution)