Skip to content

Commit

Permalink
Automatic readme update
Browse files Browse the repository at this point in the history
  • Loading branch information
actions-user authored and jamesmbaazam committed Oct 6, 2023
1 parent 3331451 commit c19947c
Show file tree
Hide file tree
Showing 2 changed files with 313 additions and 1 deletion.
314 changes: 313 additions & 1 deletion README.md
Original file line number Diff line number Diff line change
Expand Up @@ -59,7 +59,319 @@ library("epichains")

# Quick start

Work in progress
## Core functionality

*epichains* provides four main functions:

### [likelihood()](https://epiverse-trace.github.io/epichains/reference/likelihood.html)

This function calculates the likelihood/loglikelihood of observing a
vector of outbreak summaries obtained from transmission chains.
Summaries here refer to transmission chain sizes or lengths/durations.

`likelihood()` requires a vector of chain summaries (sizes or lengths),
`chains`, the corresponding statistic to calculate, `statistic`, and the
offspring distribution, `offspring_dist` its associated parameters. It
also requires `nsim_obs`, which is the number of simulations to run if
the likelihoods do not have a closed-form solution and must be
simulated. This argument will be explained further in the [“Getting
Started”](https://epiverse-trace.github.io/epichains/articles/epichains.html)
vignette.

Let’s look at the following example where we estimate the loglikelihood
of observing `chain_sizes`.

``` r
set.seed(121)
# example of observed chain sizes
# randomly generate 20 chains of size between 1 to 10
chain_sizes <- sample(1:10, 20, replace = TRUE)
```

``` r
# estimate loglikelihood of the observed chain sizes
likelihood_eg <- likelihood(
chains = chain_sizes,
statistic = "size",
offspring_dist = "pois",
nsim_obs = 100,
lambda = 0.5
)
# Print the estimate
likelihood_eg
#> [1] -67.82879
```

### [simulate_tree()](https://epiverse-trace.github.io/epichains/reference/simulate_tree.html)

`simulate_tree()` simulates an outbreak from a given number of
infections. It retains and returns information on infectors (ancestors),
infectees, the generation of infection, and the time, if a serial
distribution is specified.

Let’s look at an example where we simulate the transmission trees of
$10$ initial infections/chains. We assume a poisson offspring
distribution with mean, $\text{lambda} = 0.9$, and a serial interval of
$3$ days:

``` r
set.seed(123)

sim_tree_eg <- simulate_tree(
nchains = 10,
statistic = "size",
offspring_dist = "pois",
stat_max = 10,
serials_dist = function(x) 3,
lambda = 0.9
)

head(sim_tree_eg)
#> < tree head (from first known ancestor) >
#> chain_id sim_id ancestor generation time
#> 11 2 2 1 2 3
#> 13 3 2 1 2 3
#> 14 4 2 1 2 3
#> 16 5 2 1 2 3
#> 19 7 2 1 2 3
#> 20 8 2 1 2 3
```

`simulate_tree()` can model population-level intervention by reducing
the $R_0$, using the `intvn_mean_reduction` argument.

To illustrate this, we will use the previous example and specify a
population-level intervention that reduces $R_0$ by $50\%$.

``` r
set.seed(123)

sim_tree_intvn_eg <- simulate_tree(
nchains = 10,
statistic = "size",
offspring_dist = "pois",
intvn_mean_reduction = 0.5,
stat_max = 10,
serials_dist = function(x) 3,
lambda = 0.9
)

head(sim_tree_intvn_eg)
#> < tree head (from first known ancestor) >
#> chain_id sim_id ancestor generation time
#> 11 2 2 1 2 3
#> 12 4 2 1 2 3
#> 13 5 2 1 2 3
#> 15 8 2 1 2 3
#> 14 5 3 1 2 3
#> 16 2 3 2 3 6
```

### [simulate_summary()](https://epiverse-trace.github.io/epichains/reference/simulate_summary.html)

`simulate_summary()` is basically `simulate_tree()` except that it does
not retain information on each infector and infectee. It returns the
eventual size or length/duration of each transmission chain.

Here is an example to simulate the previous examples without
intervention, returning the size of each of the $10$ chains. It assumes
a poisson offspring distribution with mean of $0.9$.

``` r
set.seed(123)

simulate_summary_eg <- simulate_summary(
nchains = 10,
statistic = "size",
offspring_dist = "pois",
stat_max = 10,
lambda = 0.9
)

# Print the results
simulate_summary_eg
#> `epichains` object
#>
#> [1] 1 Inf 4 4 Inf 1 2 Inf 5 3
#>
#> Simulated chain sizes:
#>
#> Max: 5
#> Min: 1
```

Here is an example with an intervention that reduces $R_0$ by $50\%$.

``` r
simulate_summary_intvn_eg <- simulate_summary(
nchains = 10,
statistic = "size",
offspring_dist = "pois",
intvn_mean_reduction = 0.5,
stat_max = 10,
lambda = 0.9
)

# Print the results
simulate_summary_intvn_eg
#> `epichains` object
#>
#> [1] 1 1 1 1 1 2 1 1 2 1
#>
#> Simulated chain sizes:
#>
#> Max: 2
#> Min: 1
```

### [simulate_tree_from_pop()](https://epiverse-trace.github.io/epichains/reference/simulate_tree_from_pop.html)

`simulate_tree_from_pop()` simulates outbreaks based on a specified
population size and pre-existing immunity until the susceptible pool
runs out.

Here is a quick example where we simulate an outbreak in a population of
size $1000$. We assume individuals have a poisson offspring distribution
with mean, $\text{lambda} = 1$, and serial interval of $3$:

``` r
set.seed(7)

sim_tree_from_pop_eg <- simulate_tree_from_pop(
pop = 1000,
offspring_dist = "pois",
lambda = 1,
serials_dist = function(x) {3}
)

head(sim_tree_from_pop_eg)
#> < tree head (from first known ancestor) >
#> sim_id ancestor generation time
#> 2 2 1 2 3
#> 3 3 1 2 3
#> 4 4 1 2 3
#> 5 5 1 2 3
#> 6 6 2 3 6
#> 7 7 6 4 9
```

## Other functionalities

### Summarising

You can run `summary()` on `<epichains>` objects to get useful
summaries.

``` r
# Example with simulate_tree()
set.seed(123)

sim_tree_eg <- simulate_tree(
nchains = 10,
statistic = "size",
offspring_dist = "pois",
stat_max = 10,
serials_dist = function(x) 3,
lambda = 0.9
)

summary(sim_tree_eg)
#> $chains_ran
#> [1] 10
#>
#> $max_time
#> [1] 12
#>
#> $unique_ancestors
#> [1] 9
#>
#> $max_generation
#> [1] 5

# Example with simulate_summary()
set.seed(123)

simulate_summary_eg <- simulate_summary(
nchains = 10,
statistic = "size",
offspring_dist = "pois",
stat_max = 10,
lambda = 0.9
)

# Get summaries
summary(simulate_summary_eg)
#> $chain_ran
#> [1] 10
#>
#> $max_chain_stat
#> [1] 5
#>
#> $min_chain_stat
#> [1] 1
```

### Aggregating

You can aggregate `<epichains>` objects returned by the `simulate_*()`
functions into a time series, which is a `<data.frame>` with columns
“cases” and either “generation” or “time”, depending on the value of
`grouping_var`.

To aggregate over “time”, you must have specified a serial interval
distribution in the simulation step.

``` r
# Example with simulate_tree()
set.seed(123)

sim_tree_eg <- simulate_tree(
nchains = 10,
statistic = "size",
offspring_dist = "pois",
stat_max = 10,
serials_dist = function(x) 3,
lambda = 0.9
)

aggregate(sim_tree_eg, grouping_var = "time")
#> time cases
#> 1 0 10
#> 2 3 13
#> 3 6 15
#> 4 9 18
#> 5 12 2
```

### Plotting

Aggregated `<epichains>` objects can easily be plotted using base R or
`ggplot2` with little to no data manipulation.

Here is an end-to-end example from simulation through aggregation to
plotting.

``` r
# Run simulation with simulate_tree()
set.seed(123)

sim_tree_eg <- simulate_tree(
nchains = 10,
statistic = "size",
offspring_dist = "pois",
stat_max = 10,
serials_dist = function(x) 3,
lambda = 0.9
)

# Aggregate cases over time
sim_aggreg <- aggregate(sim_tree_eg, grouping_var = "time")

# Plot cases over time
plot(sim_aggreg, type = "b")
```

<img src="man/figures/README-unnamed-chunk-13-1.png" width="100%" />

## Package vignettes

Expand Down
Binary file added man/figures/README-unnamed-chunk-13-1.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.

0 comments on commit c19947c

Please sign in to comment.