Full Workflow Tutorial

In this example, we will use a baseline simulated dataset to show a full workflow using the EpiFusion Framework: (i) data preparation (ii) prior and parameter specification (iii) running EpiFusion (iv) parsing and plotting output. This tutorial is also outlined in the F1000 Publication 'The EpiFusion Analysis Framework for joint phylodynamic and epidemiological analysis of outbreak characteristics'

We will use the EpiFusionUtilities R package, which can be installed using devtools:

devtools::install_github("https://github.com/ciarajudge/EpiFusionUtilities")

Dataset

Data Preparation

First we load and inspect the data for this example using the EpiFusionUtilities function baseline_dataset(). This function loads a data frame with weekly case incidence (formatted with two columns, Cases and Date), a fixed time-scaled phylogenetic tree of samples, and samples from a tree posterior (with burn-in removed) from a BEAST analysis which we will use in a different tutorial (Phylogenetic Uncertainty).

baseline_dataset()

Next we set two date objects: the 'index date', or the earliest date from which we will model the outbreak origin date, and the date of sampling of the last observed sequence from the dataset. Whilst for this example we know (through the simulation process) that the outbreak origin was the 1st of January 2024, it can be good practice to set the index date to some time before the date that we suspect the outbreak began in the location represented by our case and phylogenetic data, to ensure the dynamics of full outbreak are captured.

index_date <- as.Date("2023-12-26")
last_sequence <- as.Date("2024-03-10")

To prepare the tree objects for EpiFusion we can use the prepare_epifusion_tree function from EpiFusionUtilities. This function processes the tree(s) for input to EpiFusion and writes them to the provided file path. In the case where a single fixed tree (our first example) is provided to this function it also returns the processed tree as an R phylo object, which here we reassign to the variable fixed_tree.

fixed_tree <- prepare_epifusion_tree(baseline_tree, index_date, last_sequence, "Data/Processed/baseline_fixed_tree.tree")

Definition of parameters

We will create an EpiFusion XML file using the generate_epifusion_xml function from EpiFusionUtilities. This function populates the below XML template with our data and creates a new file, and other arguments to this function can adjust parameters from their default value.

<?xml version="1.0" encoding="UTF-8"?>
<EpiFusionInputs>
  <loggers>
    <fileBase>FILESTEM</fileBase>
    <logEvery>10</logEvery>
  </loggers>
  <data>
    <incidence>
      <incidenceVals>INCIDENCE</incidenceVals>
      <incidenceTimes type="exact">INCIDENCETIMES</incidenceTimes>
    </incidence>
    <tree>
      <treePosterior></treePosterior>
    </tree>
    <epicontrib>0.5</epicontrib>
    <changetimes>0</changetimes>
  </data>
  <analysis>
    <type>looseformbeta</type>
    <startTime>null</startTime>
    <endTime>null</endTime>
    <inferTimeOfIntroduction>false</inferTimeOfIntroduction>
  </analysis>
  <model>
    <epiObservationModel>poisson</epiObservationModel>
  </model>
  <parameters>
    <epiOnly>false</epiOnly>
    <phyloOnly>false</phyloOnly>
    <numParticles>200</numParticles>
    <numSteps>2000</numSteps>
    <numThreads>8</numThreads>
    <numChains>4</numChains>
    <stepCoefficient>0.05</stepCoefficient>
    <resampleEvery>7</resampleEvery>
    <segmentedDays>true</segmentedDays>
    <samplingsAsRemovals>1</samplingsAsRemovals>
    <pairedPsi>false</pairedPsi>
  </parameters>
  <priors>
    <gamma>
      <stepchange>false</stepchange>
      <disttype>TruncatedNormal</disttype>
      <mean>0.143</mean>
      <standarddev>0.05</standarddev>
      <lowerbound>0.0</lowerbound>
    </gamma>
    <psi>
      <stepchange>false</stepchange>
      <disttype>TruncatedNormal</disttype>
      <mean>0.001</mean>
      <standarddev>0.0005</standarddev>
      <lowerbound>0.0</lowerbound>
    </psi>
    <phi>
      <stepchange>false</stepchange>
      <disttype>TruncatedNormal</disttype>
      <mean>0.02</mean>
      <standarddev>0.01</standarddev>
      <lowerbound>0.0</lowerbound>
    </phi>
    <initialBeta>
      <stepchange>false</stepchange>
      <disttype>Uniform</disttype>
      <min>0.3</min>
      <max>0.8</max>
    </initialBeta>
    <betaJitter>
      <stepchange>false</stepchange>
      <disttype>Uniform</disttype>
      <min>0.001</min>
      <max>0.05</max>
    </betaJitter>
  </priors>
</EpiFusionInputs>

We will generate an EpiFusion XML using the fixed tree we prepared with the prepare_epifusion_tree function and our loaded case incidence data. First we will make lists of the various parts of the XML file we wish to override from the default. For example, the below code represents the loggers chunk in the default XML that details how often we sample from the MCMC (every 10 MCMC steps):

  <loggers>
    <fileBase>FILESTEM</fileBase>
    <logEvery>10</logEvery>
  </loggers>

To override this, we will make a list in R that we will later pass to the loggers argument of the generate_epifusion_xml function to specify our output folder filepath as Results/fixed_tree and sample from the MCMC chain every 5 steps:

loggers <- list(fileBase = "Results/baseline_fixed_tree", logEvery = 5)

We will also slightly adjust the prior for initialBeta, or $\beta_0$ (force of infection at the beginning of the time series) from the default settings. As the default prior for $\gamma$ is a truncated normal distribution with mean $0.15$, by setting the initial $\beta$ value as $0.1 < \beta < 0.5$ we indicate that the initial $R_t$ is approximately between $0.66$ and $3.33$.

priors <- list(initialBeta = list(stepchange = "false",
                            disttype = "Uniform",
                            min = 0.1,
                            max = 0.5))

In this example we are happy with the other parameters in the default XML, so we can generate the XML file Data/EpiFusion_XMLs/fixed_tree_inputfile.xml with the following code:

generate_epifusion_XML(tree = "Data/Processed/baseline_fixed_tree.tree",
                       case_incidence = baseline_caseincidence,
                       index_date = index_date,
                       loggers = loggers,
                       priors = priors,
                       xml_filepath = "Data/EpiFusion_XMLs/baseline_fixed_tree_inputfile.xml")

Running EpiFusion

To run EpiFusion for the fixed tree example, we will use the run_epifusion function from EpiFusionUtilities to run the program within our R session:

run_epifusion("Data/EpiFusion_XMLs/baseline_fixed_tree_inputfile.xml")

On conclusion of its analysis, EpiFusion saves a timings.txt file to the output folder with the total runtime in nanoseconds, which we examine and convert to minutes below:

runtime <- suppressWarnings(read.table("Results/baseline_fixed_tree/timings.txt")[1,1]) / 6e10
paste0("Runtime: ",runtime," minutes")

Parsing and plotting the output

First we will use the load_raw_epifusion function to import the full raw results. This function automatically produces plots (Figure 4) of the likelihood and parameter traces using the plot_likelihood_trace and plot_parameter_trace functions (these plots can be suppressed by passing the argument suppress = TRUE to the function). This allows us to check for convergence and help to identify what proportion of each chain to discard as burn-in.

raw_output_fixed <- load_raw_epifusion("Results/baseline_fixed_tree/")

Next we can discard the burn-in from each MCMC chain and combine all chains into a combined posterior using the extract_posterior_epifusion function which takes a raw EpiFusion object and the proportion of each chain to discard as burn-in as its arguments. By default, the function returns means and Highest Posterior Density (HPD) intervals for the trajectories and parameters fitted by EpiFusion, however by specifying include_samples = TRUE we also instruct the function to return the actual posterior samples (minus burn-in) for inspection. This greatly increases the memory used by the posterior output object in your R environment, so is recommended for initial inspection of your results but not for downstream tasks such as loading posteriors from many analyses for plotting.

parsed_output_fixed <- extract_posterior_epifusion(raw_output_fixed, 0.1, include_samples = TRUE)

The extracted posterior object from the extract_posterior_epifusion function contains mean and HPD intervals of increasing width for infection, R_t, cumulative infection and fitted epidemiological case trajectories. The trajectory_table function can parse these into a convenient table structured to be suitable for plotting with ggplot2.

traj_table <- trajectory_table(parsed_output_fixed, index_date)

It is possible use this table with ggplot functions to plot and inspect the inferred trajectories. However we also provide a function, plot_trajectories that takes the trajectory table as input and automatically plots all three trajectory types.

plot_trajectories(traj_table)

The plot_trajectories function also takes additional arguments to allow more customisation. For example, it is possible to provide a specific trajectory type to plot using the type argument, and specify bespoke plot colours using the plot_colours argument. Here we will plot only the $R_t$ trajectories, in pink .

plot_trajectories(traj_table, type = "rt", plot_colours = "pink")

As this was a combined analysis that has used case incidence data, it is possible to examine the fit of the case incidence simulated within the model to the provided data. We already have the case incidence data loaded from the data preparation stage, so we can add the mean and HPD intervals of the fit to the existing table.

epi_data_and_fit_table <- baseline_caseincidence %>%
  mutate(Mean_Case_Fit = parsed_output_fixed$fitted_epi_cases$mean_fitted_epi_cases,
         Lower95_Cases = parsed_output_fixed$fitted_epi_cases$fitted_epi_cases_hpdintervals$HPD0.95$Lower,
         Upper95_Cases = parsed_output_fixed$fitted_epi_cases$fitted_epi_cases_hpdintervals$HPD0.95$Upper,
         Lower88_Cases = parsed_output_fixed$fitted_epi_cases$fitted_epi_cases_hpdintervals$HPD0.88$Lower,
         Upper88_Cases = parsed_output_fixed$fitted_epi_cases$fitted_epi_cases_hpdintervals$HPD0.88$Upper,
         Lower66_Cases = parsed_output_fixed$fitted_epi_cases$fitted_epi_cases_hpdintervals$HPD0.66$Lower,
         Upper66_Cases = parsed_output_fixed$fitted_epi_cases$fitted_epi_cases_hpdintervals$HPD0.66$Upper)


ggplot(epi_data_and_fit_table, aes(x = Date)) +
  geom_line(aes(y = Mean_Case_Fit), col = "#e95b0d") +
  geom_ribbon(aes(ymin = Lower95_Cases, ymax = Upper95_Cases), fill = "#e95b0d", alpha = 0.2) +
  geom_ribbon(aes(ymin = Lower88_Cases, ymax = Upper88_Cases), fill = "#e95b0d",alpha = 0.2) +
  geom_ribbon(aes(ymin = Lower66_Cases, ymax = Upper66_Cases), fill = "#e95b0d", alpha = 0.2) +
  geom_point(aes(y = Cases), col = "#01454f") +
  lshtm_theme()

Finally we can examine the posteriors of the MCMC parameters. The posterior extraction process uses the R package stable.GR to perform gelman-rubin convergence tests on each parameter, and estimate the effective sample sizes of each. If the gelman-rubin statistic is less than 1.015 this indicates MCMC convergence. For this example we used a minimal number of chains and steps (4 chains with 2000 MCMC steps) for time efficiency, but if the MCMC has not converged it may be necessary to run each chain for longer.

print(parsed_output_fixed$parameters$gamma$rhat)
print(parsed_output_fixed$parameters$gamma$ess)

We can also view the posterior density of a parameter by plotting the samples from the MCMC, which we can access from the posterior object due to setting include_samples = TRUE when we extracted the posterior earlier using extract_epifusion_posterior.

ggplot(data = data.frame(Gamma = parsed_output_fixed$parameters$gamma$samples), aes(x = Gamma)) +
  geom_density(fill = "#01454f", alpha = 0.3) +
  lshtm_theme()