diff --git a/.Rbuildignore b/.Rbuildignore index 70503e2..e6a88c2 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -5,3 +5,8 @@ LICENSE.md codecov.yml ^doc$ ^Meta$ +^_pkgdown\.yml$ +^docs$ +^pkgdown$ +^\.github$ +^README\.Rmd$ diff --git a/.github/workflows/pkgdown.yaml b/.github/workflows/pkgdown.yaml new file mode 100644 index 0000000..bfc9f4d --- /dev/null +++ b/.github/workflows/pkgdown.yaml @@ -0,0 +1,49 @@ +# Workflow derived from https://github.com/r-lib/actions/tree/v2/examples +# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help +on: + push: + branches: [main, master] + pull_request: + release: + types: [published] + workflow_dispatch: + +name: pkgdown.yaml + +permissions: read-all + +jobs: + pkgdown: + runs-on: ubuntu-latest + # Only restrict concurrency for non-PR jobs + concurrency: + group: pkgdown-${{ github.event_name != 'pull_request' || github.run_id }} + env: + GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} + permissions: + contents: write + steps: + - uses: actions/checkout@v4 + + - uses: r-lib/actions/setup-pandoc@v2 + + - uses: r-lib/actions/setup-r@v2 + with: + use-public-rspm: true + + - uses: r-lib/actions/setup-r-dependencies@v2 + with: + extra-packages: any::pkgdown, local::. + needs: website + + - name: Build site + run: pkgdown::build_site_github_pages(new_process = FALSE, install = FALSE) + shell: Rscript {0} + + - name: Deploy to GitHub pages 🚀 + if: github.event_name != 'pull_request' + uses: JamesIves/github-pages-deploy-action@v4.5.0 + with: + clean: false + branch: gh-pages + folder: docs diff --git a/.gitignore b/.gitignore index e58875e..5d394ef 100644 --- a/.gitignore +++ b/.gitignore @@ -75,3 +75,4 @@ data-raw/ # Output from vignettes figure/ +docs diff --git a/DESCRIPTION b/DESCRIPTION index 2b01ffa..e43cc8f 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -49,3 +49,5 @@ Suggests: VignetteBuilder: knitr Config/testthat/edition: 3 LazyData: true +URL: https://traitecoevo.github.io/hmde/ +BugReports: https://github.com/traitecoevo/austraits/issues diff --git a/README.Rmd b/README.Rmd new file mode 100644 index 0000000..70ecf08 --- /dev/null +++ b/README.Rmd @@ -0,0 +1,86 @@ +--- +output: github_document +--- + +# hmde + +The goal of hmde is to fit a model for the rate of change in some quantity based on a set of pre-defined functions arising from ecological applications. We estimate differential equation parameters from repeated observations of a process, such as growth rate parameters to data of sizes over time. +In other language, `hmde` implements 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. + +A pre-print paper is available on [bioarXiv](https://doi.org/10.1101/2025.01.15.633280), or as the [hmde_paper.pdf](https://github.com/traitecoevo/hmde/blob/master/inst/doc/hmde_paper.pdf) file here. + +## 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 + +`hmde` 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. + +### 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 + +‘hmde’ is under active development. You can install the current +developmental version of ‘hmde’ from [GitHub](https://github.com/) with: + +``` r +# install.packages("remotes") +remotes::install_github("traitecoevo/hmde") +``` + +## Quick demo + +Create constant growth data with measurement error: + +``` r +library(hmde) +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 +constant_fit <- hmde_model("constant_single_ind") |> + hmde_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 + ) |> + hmde_run(chains = 1, iter = 1000, verbose = FALSE, show_messages = FALSE) +``` + +## Found a bug? + +Please submit a [GitHub issue](https://github.com/traitecoevo/hmde/issues) with details of the bug. A [reprex](https://reprex.tidyverse.org/) would be particularly helpful with the bug-proofing process! diff --git a/README.md b/README.md index daed970..746b193 100644 --- a/README.md +++ b/README.md @@ -1,68 +1,75 @@ - - -# hmde - - - -[![R-CMD-check](https://github.com/traitecoevo/hmde/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/traitecoevo/hmde/actions/workflows/R-CMD-check.yaml) -[![Codecov test -coverage](https://codecov.io/gh/traitecoevo/hmde/branch/master/graph/badge.svg)](https://app.codecov.io/gh/traitecoevo/hmde?branch=master) -[![Lifecycle: -experimental](https://img.shields.io/badge/lifecycle-experimental-orange.svg)](https://lifecycle.r-lib.org/articles/stages.html#experimental) -[![DOI: -10.1101/2025.01.15.633280](https://img.shields.io/badge/DOI%3A-10.1101%2F2025.01.15.633280-blue?link=https%3A%2F%2Fdoi.org%2F10.1101%2F2025.01.15.633280)](https://doi.org/10.1101/2025.01.15.633280) - - - -The goal of hmde is to fit a model for the rate of change in some quantity based on a set of pre-defined functions arising from ecological applications. We estimate differential equation parameters from repeated observations of a process, such as growth rate parameters to data of sizes over time. -In other language, `hmde` implements 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. - -A pre-print paper is available on [bioarXiv](https://doi.org/10.1101/2025.01.15.633280), or as the [hmde_paper.pdf](https://github.com/traitecoevo/hmde/blob/master/inst/doc/hmde_paper.pdf) file here. +# hmde + +The goal of hmde is to fit a model for the rate of change in some +quantity based on a set of pre-defined functions arising from ecological +applications. We estimate differential equation parameters from repeated +observations of a process, such as growth rate parameters to data of +sizes over time. In other language, `hmde` implements 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. + +A pre-print paper is available on +[bioarXiv](https://doi.org/10.1101/2025.01.15.633280), or as the +[hmde_paper.pdf](https://github.com/traitecoevo/hmde/blob/master/inst/doc/hmde_paper.pdf) +file here. ## 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}$$ +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. +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. +`hmde` 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. +$$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. ### 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. +$$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)) +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. - -## +$$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 @@ -79,12 +86,13 @@ remotes::install_github("traitecoevo/hmde") Create constant growth data with measurement error: ``` r +library(hmde) 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. +$$s_{ij} \sim \mathcal{N}\left( 0, \sigma_e \right)$$ blows up as +$\sigma_e$ approaches 0. Fit the model: @@ -101,4 +109,7 @@ constant_fit <- hmde_model("constant_single_ind") |> ## Found a bug? -Please submit a [GitHub issue](https://github.com/traitecoevo/hmde/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/hmde/issues) with details of the +bug. A [reprex](https://reprex.tidyverse.org/) would be particularly +helpful with the bug-proofing process! diff --git a/Read-and-delete-me b/Read-and-delete-me deleted file mode 100644 index 7ae89b5..0000000 --- a/Read-and-delete-me +++ /dev/null @@ -1,10 +0,0 @@ -Stan-specific notes: - -* All '.stan' files containing stanmodel definitions must be placed in 'inst/stan'. -* Additional files to be included by stanmodel definition files - (via e.g., #include "mylib.stan") must be placed in any subfolder of 'inst/stan'. -* Additional C++ files needed by any '.stan' file must be placed in 'inst/include', - and can only interact with the Stan C++ library via '#include' directives - placed in the file 'inst/include/stan_meta_header.hpp'. -* The precompiled stanmodel objects will appear in a named list called 'stanmodels', - and you can call them with e.g., 'rstan::sampling(stanmodels$foo, ...)' diff --git a/_pkgdown.yml b/_pkgdown.yml new file mode 100644 index 0000000..9a029e1 --- /dev/null +++ b/_pkgdown.yml @@ -0,0 +1,26 @@ +url: https://traitecoevo.github.io/hmde/ +template: + bootstrap: 5 + bootswatch: flatly + +navbar: + structure: + left: [intro, reference, articles] + right: [search, github] + components: + articles: + text: Articles + menu: + - text: "Case Studies" + - text: Constant growth model + href: 'articles/constant-growth.html' + - text: von Bertalanffy model + href: 'articles/von-bertalanffy.html' + - text: Canham model + href: 'articles/canham.html' + github: + icon: fab fa-github fa-lg + href: https://github.com/traitecoevo/hmde/ + + + diff --git a/hmde.Rproj b/hmde.Rproj index fd50596..6a3ede2 100644 --- a/hmde.Rproj +++ b/hmde.Rproj @@ -1,5 +1,4 @@ Version: 1.0 -ProjectId: 276917f4-3da1-492e-bb45-7ffd205ea048 RestoreWorkspace: No SaveWorkspace: No diff --git a/inst/figures/hex.R b/inst/figures/hex.R new file mode 100644 index 0000000..6086571 --- /dev/null +++ b/inst/figures/hex.R @@ -0,0 +1,39 @@ +library(hexSticker) + +measurement_data_transformed <- Tree_Size_Ests$measurement_data %>% + group_by(ind_id) %>% + mutate( + delta_y_obs = y_obs - lag(y_obs), + obs_interval = time - lag(time), + obs_growth_rate = delta_y_obs/obs_interval, + delta_y_est = y_hat - lag(y_hat), + est_growth_rate = delta_y_est/obs_interval + ) %>% + ungroup() + +#Plots of size over time for a sample of 5 individuals +sample_ids <- sample(1:nrow(Tree_Size_Ests$individual_data), size=5) +plot_data <- measurement_data_transformed %>% + filter(ind_id %in% sample_ids) + +plot <- ggplot(data=plot_data, aes(group = ind_id)) + + geom_point(aes(x = time, y=y_obs, colour = as.factor(ind_id)), + shape = 1) + + geom_line(aes(x = time, y=y_obs, colour = as.factor(ind_id)), + linetype = "dashed") + + geom_point(aes(x = time, y=y_hat, colour = as.factor(ind_id)), + shape = 2) + + geom_line(aes(x = time, y=y_hat, colour = as.factor(ind_id)), + linetype = "solid") + + labs(x="Time (years)", y="DBH (cm)", colour="Ind. ID") + + theme_void() + + theme(legend.position = "none") + +hexSticker::sticker(plot, package="hmde", + p_color = "white", + p_size=40, + p_y = 1, + s_x=1, s_y=1, + s_width=1.35, s_height = 1.35, + h_fill = "navy", h_color = "white", + filename="inst/figures/hmde_hex.png") diff --git a/inst/figures/hmde_hex.png b/inst/figures/hmde_hex.png new file mode 100644 index 0000000..7abf3bc Binary files /dev/null and b/inst/figures/hmde_hex.png differ diff --git a/man/figures/hex.R b/man/figures/hex.R new file mode 100644 index 0000000..6086571 --- /dev/null +++ b/man/figures/hex.R @@ -0,0 +1,39 @@ +library(hexSticker) + +measurement_data_transformed <- Tree_Size_Ests$measurement_data %>% + group_by(ind_id) %>% + mutate( + delta_y_obs = y_obs - lag(y_obs), + obs_interval = time - lag(time), + obs_growth_rate = delta_y_obs/obs_interval, + delta_y_est = y_hat - lag(y_hat), + est_growth_rate = delta_y_est/obs_interval + ) %>% + ungroup() + +#Plots of size over time for a sample of 5 individuals +sample_ids <- sample(1:nrow(Tree_Size_Ests$individual_data), size=5) +plot_data <- measurement_data_transformed %>% + filter(ind_id %in% sample_ids) + +plot <- ggplot(data=plot_data, aes(group = ind_id)) + + geom_point(aes(x = time, y=y_obs, colour = as.factor(ind_id)), + shape = 1) + + geom_line(aes(x = time, y=y_obs, colour = as.factor(ind_id)), + linetype = "dashed") + + geom_point(aes(x = time, y=y_hat, colour = as.factor(ind_id)), + shape = 2) + + geom_line(aes(x = time, y=y_hat, colour = as.factor(ind_id)), + linetype = "solid") + + labs(x="Time (years)", y="DBH (cm)", colour="Ind. ID") + + theme_void() + + theme(legend.position = "none") + +hexSticker::sticker(plot, package="hmde", + p_color = "white", + p_size=40, + p_y = 1, + s_x=1, s_y=1, + s_width=1.35, s_height = 1.35, + h_fill = "navy", h_color = "white", + filename="inst/figures/hmde_hex.png") diff --git a/man/figures/hmde_hex.png b/man/figures/hmde_hex.png new file mode 100644 index 0000000..7abf3bc Binary files /dev/null and b/man/figures/hmde_hex.png differ diff --git a/man/hmde-package.Rd b/man/hmde-package.Rd index 7e86a43..90e7348 100644 --- a/man/hmde-package.Rd +++ b/man/hmde-package.Rd @@ -10,6 +10,14 @@ A package to implement a selection of hierarchical Bayesian longitudinal models } \references{ Stan Development Team (NA). RStan: the R interface to Stan. R package version 2.26.23. https://mc-stan.org +} +\seealso{ +Useful links: +\itemize{ + \item \url{https://traitecoevo.github.io/hmde/} + \item Report bugs at \url{https://github.com/traitecoevo/austraits/issues} +} + } \author{ \strong{Maintainer}: Tess O'Brien \email{theresa.obrien@unsw.edu.au} (\href{https://orcid.org/XXXX-XXXX-XXXX-XXXX}{ORCID}) [copyright holder]