Extract underlying distributions and plot to compare

[adamkucharski]

From discussions, thought it would be useful for user have ability to plot the available distributions and have these stored for future use as epidist objects. Below code extracts onset-to-death parameters for Ebola (so opening steps similar to quick start in this discusison ), then converts distributions where available and plots and stores.

The resulting epidist_ebola_otd_list object could then be used as model input distribution, e.g. in simulation scenarios like the one described here.

library(epireview)
library(epiparameter)

# Load Ebola data
data_in <- load_epidata(pathogen = "ebola")

# Extract parameters for onset to death, with qualtiy score > 50
params_in <- data_in[["params"]] |> dplyr::filter(article_qa_score>50)
params_si <- params_in |> dplyr::filter(parameter_type=="Human delay - Symptom Onset/Fever to Death") #"Human delay - serial interval")

# Select distribution
params_si_with_dist <- params_si |> dplyr::filter(!is.na(distribution_type))

# Start by extracting only those with Gamma
# NOTE: Need a check here to iterate across all distributions
params_si_with_dist <- params_si |> dplyr::filter(distribution_type=="Gamma")

# Create epidist objects for each one
epidist_ebola_otd_list <- list()

par(mar=c(3,3,1,1),mgp=c(2,0.7,0))

x_range <- seq(0,20,0.1)
plot(1, type="n", xlab="onset to death (days)", ylab="pdf", xlim=c(0,max(x_range)), ylim=c(0,0.12))

for(ii in 1:nrow(params_si_with_dist)){
  
  # NOTE: need check on parameter types here - currently assuming mean/sd
  # Can use `params_si_with_dist$distribution_par1_type`
  
  mean_ii <- as.numeric(params_si_with_dist$distribution_par1_value[ii])
  sd_ii <- as.numeric(params_si_with_dist$distribution_par2_value[ii])
  score_ii <- as.numeric(params_si_with_dist$article_qa_score[ii])
  
  # Get parameters
  dist_params <- epiparameter::convert_summary_stats_to_params(
    "gamma",
    mean = mean_ii,
    sd = sd_ii
  )
  
  # Create epidist object
  epidist_ebola_otd = epiparameter::epidist(
    "Ebola virus disease",
    "Ebolavirus",
    epi_dist = "onset_to_death",
    prob_distribution = "gamma",
    prob_distribution_params = unlist(dist_params)
  )
  
  # Store output in list
  epidist_ebola_otd_list[[ii]] <- epidist_ebola_otd
  
  # Plot distribution, with colour intensity represented by score
  lines(x_range,density(epidist_ebola_otd,x_range),lwd=2,col=rgb(0,0.5,1,score_ii/100))

}

[adamkucharski]

Collaborative design document for function(s) to extract ‘distribution of distributions’ from epireview using epiparameter

1. Extract available studies and convert to distributions
Several different possibilities we would have to deal with to generate a distribution of distributions (from easier to harder):
a) Known distribution (e.g. Gamma) and available parameters (e.g. distribution_par1_type and distribution_par2_type). As outlined in discussion, can use epiparameter to convert to distribution parameter (e.g. with convert_summary_stats_to_params) and store as epidist object (e.g. list?).

b) No known distribution but summary parameters available (e.g. parameter_value_type and parameter_uncertainty_type); can be linked by study id. Would require distribution assumption to convert using epiparameter (e.g. extract_param(type = "percentiles",values = c(1, 10),distribution = "gamma",percentiles = c(0.025, 0.975))). Perhaps use the most populous distribution?
Could also use the $models component of the dataset to identify estimates derived from SIR models, which would likely have be exponentially distributed?

c) Range of values with no known distribution (type specified in parameter_type) across levels of disaggregation (parameter_lower_bound and parameter_upper_bound). Can extract with epiparameter (e.g. extract_param(type = "range",values = c(10, 5, 15),distribution = "lnorm",samples = 25) ).

Additional considerations:

• How to make it easy for user to use epiparameter input (e.g. onset_to_death) and extracting relevant entries (e.g. "Human delay - Symptom Onset/Fever to Death"). Make use of emerging get_si etc. functionality in epireview?
• Need some linkage with metadata (e.g. context, sample size). Perhaps simplest to store parameter_data_id to allow user to go back to original epireview. But also value to having key values (e.g. sample size) stored in the epidist object.
• Need to ensure citation info is correctly transferred across.
• Would need option for user to decide how much of the data they want to use vs how many assumptions being made.

2. Synthesise output (i.e. list with distribution of distributions) into meaningful function for use in analysis
a) Some approaches to make use of distribution of distributions:
Estimate joint distribution (where feasible, for distributions that has analytical expression for linear combination)

b) Define a new epidist object or define a new empirical distribution (i.e. mixture model via ecdf) then (b.1) consider fitting to this (but issues with the tail potentially)
JWL: Seems very plausible to use <ecdf> which is in base R: https://stat.ethz.ch/R-manual/R-devel/library/stats/html/ecdf.html

c) Create function to generate random sample from the combined empirical distribution.
JWL: What is meant by combined distributions here? AK: same as (b), just generating samples from each distribution and sampling in weighted way rather than defining full empirical distribution.

(b) and (c) perhaps most straightforward and general, with easy intuition for users. For all approaches would need to weight mixture somehow (e.g. by sample size or inverse variance)

Broader issues to consider:

• If we want infection-to-death, may need to extract incubation period and onset-to-death. Is there easy way to do this in consistent way (rather than extracting separately and combining manually)?
• Where do functions live? To keep dependencies light, probably better to have functionality to build distribution of distributions in epiparameter, because don’t require epireview dependency (just a pre-loaded object). But useful to have the filtering options in the epireview package, especially as can then propagate through.

Original document: https://docs.google.com/document/d/18JNyoHkbIIw96rha2oEMF2DT839q-3I05fchbI2iwpU/edit?usp=sharing