Skip to content

Commit

Permalink
Update readme (#16)
Browse files Browse the repository at this point in the history
Added a general description to the README as well as some mathematical details about the included models.

* Started writing the README. Current sections describe the overall purpose and the maths structure with a short description of each inbuilt model.

* Did some further editing of the README file.

* Added links and a citation to the readme. Updated the model demo to show how to extract the estimate of constant growth rate.

* Removed .rmd file as the maths wasn't rendering properly.

* Changed line breaks to hopefully encourage the maths to render.

* Testing bracketing macro.

* Corrected bracketing.

* Had problems with '\r' turning into a new line.

* Had problems with _i being italics.

* Missed some characters.

* Missed some underscores.

* Testing white space.

* Testing white space for final equations.

* Testing \bigg vs \left.

* Correcting size and direction macros for parens.

* Removed unnecessary bit.

* Typo in code for demo.

* Modified whitespace in code demo to look nicer.

* Corrected typo for beta.

* Changed y_0_obs to be the first obs with measurement error rather than just a number.

* Update README.md

 Removed a backslash that I thought I previously got rid of.
  • Loading branch information
Tess-LaCoil authored Apr 22, 2024
1 parent e049daf commit 100074f
Show file tree
Hide file tree
Showing 2 changed files with 66 additions and 58 deletions.
48 changes: 0 additions & 48 deletions README.Rmd

This file was deleted.

76 changes: 66 additions & 10 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,51 @@ coverage](https://codecov.io/gh/traitecoevo/rmot/branch/master/graph/badge.svg)]
experimental](https://img.shields.io/badge/lifecycle-experimental-orange.svg)](https://lifecycle.r-lib.org/articles/stages.html#experimental)
<!-- badges: end -->

The goal of ‘rmot’ is to …
The goal of rmot is to implement hierarchical Bayesian longitudinal models to solve the Bayesian inverse problem of estimating differential equation parameters based on repeat measurement surveys. Estimation is done using Markov Chain Monte Carlo, implemented through
[Stan](https://mc-stan.org/) via [RStan](https://mc-stan.org/users/interfaces/rstan), built under [R](https://cran.r-project.org/) 4.3.3. The inbuilt models are based on case studies in ecology.

## The Maths

The general use case is to estimate a vector of parameters $\boldsymbol{\theta}$ for a chosen differential equation
$$f\left( Y \left( t \right), \boldsymbol{\theta} \right) = \frac{dY}{dt}$$
based on the longitudinal structure
$$Y \left( t_{j+1} \right) = Y\left( t_j \right) + \int_{t_j}^{t_{j+1}}f\left( Y \left( t \right), \boldsymbol{\theta} \right)\,dt. $$

The input data are observations of the form $y_{ij}$ for individual $i$ at time $t_j$, with repeated observations coming from the same individual. We parameterise $f$ at the individual level by estimating $\boldsymbol{\theta}\_i$ as the vector of parameters. We have hyper-parameters that determine the distribution of $\boldsymbol{\theta}\_i$ with typical prior distribution
$$\boldsymbol{\theta}\_i \sim \log \mathcal{N}\left(\boldsymbol{\mu}\_{\log\left(\boldsymbol{\theta}\right)}, \boldsymbol{\sigma}\_{\log \left( \boldsymbol{\theta} \right)}\right), $$
where $\boldsymbol{\mu}\_{\log\left(\boldsymbol{\theta}\right)}$ and $\boldsymbol{\sigma}\_{\log\left(\boldsymbol{\theta}\right)}$ are vectors of means and standard deviations. In the case of a single individual, these are chosen prior values. In the case of a multi-individual model $\boldsymbol{\mu}\_{\log\left(\boldsymbol{\theta}\right)}$ and $\boldsymbol{\sigma}\_{\log\left(\boldsymbol{\theta}\right)}$ have their own prior distributions and are fit to data.

## Implemented Models

Rmot comes with four DEs built and ready to go, each of which has a version for a single individual and multiple individuals.

### Constant Model

The constant model is given by
$$f \left( Y \left( t \right), \beta \right) = \frac{dY}{dt} = \beta,$$
and is best understood as describing the average rate of change over time.

### Power law

The power law model is given by
$$f \left( Y \left( t \right), \beta_0, \beta_1, \bar{Y} \right) = \frac{dY}{dt} = \beta_0 \bigg( \frac{Y \left( t \right)}{\bar{Y}} \bigg)^{\beta_1},$$
where $\beta_0>0$ is the coefficient, $\beta_1$ is the power, and $\bar{Y}$ is a user-provided parameter that centres the model in order to avoid correlation between the $\beta$ parameters.

### von Bertalanffy

The von Bertalanffy mode is given by
$$f \left( Y \left( t \right), \beta, Y_{max} \right) = \frac{dY}{dt} = \beta \left( Y_{max} - Y \left( t \right) \right),$$
where $\beta$ is the growth rate parameter and $Y_{max}$ is the maximum value that $Y$ takes.

### Canham

The Canham ([Canham et
al. 2004](https://doi.org/10.1890/1051-0761(2006)016%5B0540:NAOCTC%5D2.0.CO;2))
model is a hump-shaped function given by
$$f \left( Y \left( t \right), f_{max}, Y_{max}, k \right) = \frac{dY}{dt} = f_{max} \exp \Bigg( -\frac{1}{2} \bigg( \frac{ \ln \left( Y \left( t \right) / Y_{max} \right) }{k} \bigg)^2 \Bigg), $$
where $f_{max}$ is the maximum growth rate, $Y_{max}$ is the $Y$-value at which that maximum occurs, and $k$ controls how narrow or wide the peak is.

##

## Installation

Expand All @@ -26,17 +70,29 @@ remotes::install_github("traitecoevo/rmot")

## Quick demo

Create constant growth data with measurement error:

``` r
y_obs <- seq(from=2, to=15, length.out=10) + rnorm(10, 0, 0.1)
```

Measurement error is necessary as otherwise the normal likelihood
$$s_{ij} \sim \mathcal{N}\left( 0, \sigma_e \right)$$
blows up as $\sigma_e$ approaches 0.

Fit the model:

``` r
rmot_model("linear") |>
rmot_assign_data(X = Loblolly$age,
Y = Loblolly$height,
N = nrow(Loblolly)) |>
rmot_run()
constant_fit <- rmot_model("constant_single_ind") |>
rmot_assign_data(n_obs = 10, #Integer
y_obs = y_obs, #vector length n_obs
obs_index = 1:10, #vector length n_obs
time = 0:9, #Vector length n_obs
y_0_obs = y_obs[1] #Real
) |>
rmot_run(chains = 1, iter = 1000, verbose = FALSE, show_messages = FALSE)
```

## Found a bug?

Please submit a [GitHub
issue](https://github.com/traitecoevo/rmot/issues) with details of the
bug. A [reprex](https://reprex.tidyverse.org/) would be particularly
helpful with the bug-proofing process!
Please submit a [GitHub issue](https://github.com/traitecoevo/rmot/issues) with details of the bug. A [reprex](https://reprex.tidyverse.org/) would be particularly helpful with the bug-proofing process!

0 comments on commit 100074f

Please sign in to comment.