Skip to content

Commit

Permalink
Merge branch 'main' into clarify_chain_ll_docs
Browse files Browse the repository at this point in the history
  • Loading branch information
jamesmbaazam authored Jun 20, 2023
2 parents fd3abe7 + 06ed75a commit b824631
Show file tree
Hide file tree
Showing 15 changed files with 1,017 additions and 291 deletions.
1 change: 1 addition & 0 deletions .Rbuildignore
Original file line number Diff line number Diff line change
Expand Up @@ -14,3 +14,4 @@
^pkgdown$
^data-raw$
^CITATION\.cff$
^_pkgdown\.yml$
9 changes: 7 additions & 2 deletions .github/workflows/render_readme.yml
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,11 @@ jobs:
run: |
git config --local user.email "action@github.com"
git config --local user.name "GitHub Action"
git add README.md
git add README.md
# Also add README figures if they exist
if [ -d man/figures ]
then
git add man/figures/
fi
git diff-index --quiet HEAD || git commit -m "Automatic readme update"
git push origin || echo "No changes to push"
git push origin || echo "No changes to push"
1 change: 1 addition & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -31,3 +31,4 @@ rsconnect/
/Meta/
/docs/
.DS_Store
docs
46 changes: 15 additions & 31 deletions DESCRIPTION
Original file line number Diff line number Diff line change
@@ -1,42 +1,25 @@
Package: bpmodels
Title: Analysing transmission chain statistics using branching process models
Version: 0.2.1
Title: Analysing transmission chain statistics using branching process
models
Version: 0.3.0
Authors@R: c(
person(
given = "Sebastian",
family = "Funk",
email = "sebastian.funk@lshtm.ac.uk",
role = "aut",
comment = c(ORCHID = "https://orcid.org/0000-0002-2842-3406")
),
person(
given = "Zhian N.",
family = "Kamvar",
email = "zkamvar@gmail.com",
role = "ctb",
comment = c(ORCHID = "https://orcid.org/0000-0003-1458-7108")
),
person(
given = "Flavio",
family = "Finger",
email = "flavio.finger@epicentre.msf.org",
role = "aut",
comment = c(ORCHID = "https://orcid.org/0000-0002-8613-5170")
),
person(
given = "James M.",
family = "Azam",
email = "james.azam@lshtm.ac.uk",
role = c("aut", "cre"),
comment = c(ORCHID = "https://orcid.org/0000-0001-5782-7330"))
person("Sebastian", "Funk", , "sebastian.funk@lshtm.ac.uk", role = "aut",
comment = c(ORCID = "https://orcid.org/0000-0002-2842-3406")),
person("Zhian N.", "Kamvar", , "zkamvar@gmail.com", role = "ctb",
comment = c(ORCID = "https://orcid.org/0000-0003-1458-7108")),
person("Flavio", "Finger", , "flavio.finger@epicentre.msf.org", role = "aut",
comment = c(ORCID = "https://orcid.org/0000-0002-8613-5170")),
person("James M.", "Azam", , "james.azam@lshtm.ac.uk", role = c("aut", "cre"),
comment = c(ORCID = "https://orcid.org/0000-0001-5782-7330"))
)
Description: Provides methods to analyse and simulate the size and length
of branching processes with an arbitrary offspring distribution. These
can be used, for example, to analyse the distribution of chain sizes
or length of infectious disease outbreaks, as discussed in Farrington
et al. (2003) <doi:10.1093/biostatistics/4.2.279>.
License: MIT + file LICENSE
URL: https://github.com/epiverse-trace/bpmodels, https://epiverse-trace.github.io/bpmodels/
URL: https://github.com/epiverse-trace/bpmodels,
https://epiverse-trace.github.io/bpmodels/
BugReports: https://github.com/epiverse-trace/bpmodels/issues
Depends:
R (>= 3.6.0)
Expand All @@ -50,11 +33,12 @@ Suggests:
rmarkdown,
testthat,
truncdist
Config/testthat/edition: 3
VignetteBuilder:
knitr
Remotes:
github::epiverse-trace/epiparameter
Config/Needs/website:epiverse-trace/epiversetheme
Config/testthat/edition: 3
Encoding: UTF-8
LazyData: true
Roxygen: list(markdown = TRUE)
Expand Down
50 changes: 36 additions & 14 deletions NEWS.md
Original file line number Diff line number Diff line change
@@ -1,30 +1,52 @@
# bpmodels 0.3.0

## Website

* The website has been updated to use bootstrap 5. This means it has a slightly
new look.
* Exported functions under "references" have been grouped into topics.

## Vignettes

* A new vignette has been added to compile a bibliography of papers that
apply branching processes to infectious disease epidemiology.
* The "quickstart" section of the README has been moved to the "Getting started"
section of the website.
* The README now only provides a quick description of the core functions in the
package and users are referred to the website for a quick start.

## README

* The README no longer has a quickstart section as it has been moved to the
website. This is to give the README a minimalistic look.

# bpmodels 0.2.1

## Minor functionality change

* `chain_sim()` now throws a warning, instead of an error, when `tree` is set
to `FALSE` with `serial` also specified. We assume that providing a serial
interval means you want the tree of transmissions to be simulated,
so `chain_sim()` internally sets `tree = TRUE` and throws a warning explaining
what happened. This behaviour should not break any simulations with previous
versions with `bpmodels`, but if it does, please submit an issue.
To remove the warning, the user should explicitly set `tree = TRUE` when
they specify `serial`.
* `chain_sim()` now throws a warning, instead of an error, when `tree` is set
to `FALSE` with `serial` also specified. We assume that providing a serial
interval means you want the tree of transmissions to be simulated,
so `chain_sim()` internally sets `tree = TRUE` and throws a warning explaining
what happened. This behaviour should not break any simulations with previous
versions with `bpmodels`, but if it does, please submit an issue.
To remove the warning, the user should explicitly set `tree = TRUE` when
they specify `serial`.

# bpmodels 0.2.0

## Documentation

* `chain_sim()`'s help file has been updated with more details. In particular,
we describe in detail how to specify the `serial` argument as a function. We
have also added more examples.
we describe in detail how to specify the `serial` argument as a function. We
have also added more examples.

* A new vignette describing how to project COVID-19 incidence with `chain_sim()`
has been added and can be accessed on the
[bpmodels website](https://epiverse-trace.github.io/bpmodels/) under "Articles".
has been added and can be accessed on the
[bpmodels website](https://epiverse-trace.github.io/bpmodels/) under "Articles".

* The README's "quick start" section has been updated with what was
previously the introduction vignette.
* The README's "quick start" section has been updated with what was
previously the introduction vignette.

# bpmodels 0.1.9999

Expand Down
2 changes: 1 addition & 1 deletion R/borel.r
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@ dborel <- function(x, mu, log = FALSE) {
return(ld)
}

##' Generate random numbers from the Borel distribution
##' Random number generator from the Borel distribution
##'
##' Random numbers are generated by simulating from a Poisson branching process
##' @param n number of random variates to generate.
Expand Down
3 changes: 2 additions & 1 deletion R/simulate.r
Original file line number Diff line number Diff line change
Expand Up @@ -102,7 +102,8 @@
#' serial = serial_interval
#' )
#'
#' # Specifying `serial` and `tree = FALSE` will throw an error
#' # Specifying `serial` and `tree = FALSE` will throw a warning saying that
#' # `tree` was set to `TRUE` internally.
#' set.seed(123)
#' \dontrun{
#' try(chain_sim(
Expand Down
150 changes: 37 additions & 113 deletions README.Rmd
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
---
output: github_document
bibliography: vignettes/references.bib
bibliography: vignettes/references.json
link-citations: true
---

Expand Down Expand Up @@ -37,8 +37,10 @@ secondary infections caused by an infected individual.

The latest development version of the _bpmodels_ package can be installed via

```{r include=TRUE,eval=FALSE}
devtools::install_github(file.path("epiverse-trace", "bpmodels"))
```{r eval=FALSE}
# check whether {pak} is installed
if (!require("pak")) install.packages("pak")
pak::pkg_install("epiverse-trace/bpmodels")
```

To load the package, use
Expand All @@ -47,133 +49,55 @@ To load the package, use
library("bpmodels")
```

# Quick start

At the heart of the package are the `chain_ll()` and `chain_sim()` functions.

## Calculating log-likelihoods
# Core functionality

The `chain_ll()` function calculates the log-likelihood of a distribution of
chain sizes or lengths given an offspring distribution and its associated
parameters.
_bpmodels_ provides three main functions:

For example, if we have observed a distribution of chains of sizes
$1, 1, 4, 7$, we can
calculate the log-likelihood of this observed chain by assuming the offspring
per generation is Poisson distributed with a mean number (which can
be interpreted as the reproduction number $\mathcal{R_0}$) of $0.5$.

To do this, we run
`chain_ll()`: calculates the likelihoods of observing a vector of chains
of given sizes or lengths.

Here is a quick example of estimating the loglikelihood of an observed chain:
```{r}
set.seed(13)
chain_sizes <- c(1, 1, 4, 7) # example of observed chain sizes
chain_ll(x = chain_sizes, offspring = "pois", stat = "size", lambda = 0.5)
# example of observed chain sizes
chain_sizes <- c(1, 2, 3, 4)
# estimate loglikelihood of the observed chain sizes
chain_ll_eg <- chain_ll(chain_sizes, "pois", "size", lambda = 0.5)
chain_ll_eg
```

The first argument of `chain_ll()` is the chain size (or length, in number of
generations that a chain lasted) distribution to
analyse. The second argument, `offspring`, specifies the offspring
distribution. This is given as a function used to generate random offspring.
It can be any probability distribution implemented in `R`, that is, one that
has a corresponding function for generating random numbers beginning with the
letter `r`. In the case of the example above, since random Poisson numbers are
generated in `R` using a function called `rpois()`, the string to pass to the
`offspring` argument is `"pois"`.

The third argument, `stat`, determines whether to analyse chain sizes
(`"size"`, the default if this argument is not specified) or lengths
(`"length"`). Lastly, any named arguments not recognised by `chain_ll()`
are interpreted as parameters of the corresponding probability distribution,
here `lambda = 0.5` as the mean of the Poisson distribution (see the `R` help
page for the [Poisson distribution](https://stat.ethz.ch/R-manual/R-devel/library/stats/html/Poisson.html) for more information).

### Imperfect observations

By default, `chain_ll()` assumes perfect observation, where `obs_prob = 1`
(See `?chain_ll`), meaning that all transmission events are observed and
recorded in the data. If observations are imperfect, `chain_ll()` provides
the argument, `obs_prob`, for specifying the probability of observation.
This probability is used to determine the likelihood of observing the specified
chain sizes or lengths. In the case of imperfect observation, true chain sizes
or lengths are simulated repeatedly (the number of times given by the
`nsim_obs` argument), and the likelihood calculated for each of
these simulations.

For example, if the probability of observing each case is `obs_prob = 0.30`,
we use

```{r}
chain_sizes <- c(1, 1, 4, 7) # example of observed chain sizes
ll <- chain_ll(chain_sizes, "pois", "size", obs_prob = 0.3, lambda = 0.5,
nsim_obs = 10)
ll
```
`chain_sim()`: simulates transmission chains until all chains stop producing
offspring.

This returns `10` likelihood values (because `nsim_obs = 10`), which can be
averaged to come up with an overall likelihood estimate.
Below is a quick example where we simulate the chain sizes of $5$ chains with
a poisson offspring with mean, $\text{lambda} = 0.5$:
```{r}
set.seed(123)
To find out about usage of the `chain_ll()` function, you can use the `R` help
file
chain_sim_eg <- chain_sim(n = 5, offspring = "pois", stat = "size",
lambda = 0.5, tree = TRUE)
```{r eval=FALSE}
?chain_ll
head(chain_sim_eg)
```

### How `chain_ll()` works

If the probability distribution of chain sizes or lengths has an analytical
solution, this will be used. `chain_ll()` currently supports the Poisson and
negative binomial size distribution and the Poisson and geometric length
distribution.

If an analytical solution does not exist, simulations are used to approximate
this probability distributions ([using a linear approximation to the cumulative
distribution](https://en.wikipedia.org/wiki/Empirical_distribution_function)
for unobserved sizes/lengths). In that case, an extra argument `nsim_offspring`
must be passed to `chain_ll()` to specify the number of simulations to be
used for this approximation.

For example, to get offspring drawn from a binomial distribution with
probability `prob = 0.5`, we run

`chain_sim_susc()`: simulates transmission chains from a specified population
size with pre-existing immunity until the susceptible pool runs out.

Below is a quick example where we simulate chains with a poisson
offspring with mean, $\text{lambda} = 0.5$, and serial interval of $3$:
```{r}
chain_sizes <- c(1, 1, 4, 7) # example of observed chain sizes
chain_ll(chain_sizes, "binom", "size", size = 1, prob = 0.5,
nsim_offspring = 100)
```

## Simulating branching processes
set.seed(1234)
To simulate a branching process, we use the `chain_sim()` function. This
function follows the same syntax as `chain_ll()`.
chain_sim_susc_eg <- chain_sim_susc(pop = 1000, "pois",
mn_offspring = 0.5,
serial = function(x) 3
)
Below, we are simulating $5$ chains, assuming the offspring are generated using
a Poisson distribution with mean, `lambda = 0.5`. By default, `chain_sim()`
returns a vector of chain sizes/lengths. If we instead want to return
a tree of infectees and infectors, we need to specify a function for
the serial interval and set `tree = TRUE` (see next section).

```{r}
chain_sim(n = 5, offspring = "pois", stat = "size", lambda = 0.5)
head(chain_sim_susc_eg)
```

### Simulating trees

To simulate a tree of transmission chains, we specify the serial interval
generation function and set `tree = TRUE` as follows:

```{r}
set.seed(13)
serial_interval <- function(n) {
rlnorm(n, meanlog = 0.58, sdlog = 1.58)
}
chains_df <- chain_sim(
n = 5, offspring = "pois", lambda = 0.5, stat = "length",
infinite = 100, serial = serial_interval, tree = TRUE
)
head(chains_df)
```
See the ["Get started vignette"](https://epiverse-trace.github.io/bpmodels/articles/bpmodels.html) for a detailed illustration of
each function.

## Package vignettes

Expand Down
Loading

0 comments on commit b824631

Please sign in to comment.