From 095a0260e3feac4dbec8a15375dead19e5532935 Mon Sep 17 00:00:00 2001 From: cicdguy <26552821+cicdguy@users.noreply.github.com> Date: Wed, 10 May 2023 06:01:49 -0500 Subject: [PATCH 1/8] Use IDR workflows --- .github/workflows/build_docker.yaml | 67 ------------------------ .github/workflows/check.yaml | 81 +++++++++++++++++++++++++++++ .github/workflows/docs.yaml | 48 +++++++++++++++++ .github/workflows/on_push_pull.yaml | 57 -------------------- .github/workflows/release.yaml | 39 ++++++++++++++ 5 files changed, 168 insertions(+), 124 deletions(-) delete mode 100755 .github/workflows/build_docker.yaml create mode 100644 .github/workflows/check.yaml create mode 100644 .github/workflows/docs.yaml delete mode 100755 .github/workflows/on_push_pull.yaml create mode 100644 .github/workflows/release.yaml diff --git a/.github/workflows/build_docker.yaml b/.github/workflows/build_docker.yaml deleted file mode 100755 index 51bdc8b8..00000000 --- a/.github/workflows/build_docker.yaml +++ /dev/null @@ -1,67 +0,0 @@ -name: Build latest version Docker images - -env: - REGISTRY: ghcr.io - -on: - schedule: - - cron: '0 5 1 * *' - workflow_dispatch: - - -jobs: - build: - runs-on: ubuntu-latest - name: Build & Deploy Docker Images - - # Token permissions - permissions: - contents: read - packages: write - - # Build steps - steps: - - name: Checkout repository - uses: actions/checkout@v2 - - - name: Set up Docker Buildx - uses: docker/setup-buildx-action@v1 - id: buildx - with: - install: true - - - name: Log in to the Container registry - uses: docker/login-action@v1 - with: - registry: ${{ env.REGISTRY }} - username: ${{ github.actor }} - password: ${{ secrets.GITHUB_TOKEN }} - - - name: Cache Docker layers - uses: actions/cache@v2 - with: - path: /tmp/.buildx-cache - key: docker-cache-jmpost-latest - - - - name: debug - run: | - echo push = ${{ github.event_name == 'push' }} - echo tag = ${{ env.REGISTRY }}/genentech/jmpost:latest - - name: Build and push image - uses: docker/build-push-action@v2 - with: - context: ./misc/docker - push: true - tags: ${{ env.REGISTRY }}/genentech/jmpost:latest - cache-from: type=local,src=/tmp/.buildx-cache - cache-to: type=local,dest=/tmp/.buildx-cache-new - - # Fix bug where docker cache just continiously grows as it doesn't purge older layers - # That are no longer in use - # https://github.com/docker/build-push-action/blob/master/docs/advanced/cache.md#local-cache - - name: Move cache - run: | - rm -rf /tmp/.buildx-cache - mv /tmp/.buildx-cache-new /tmp/.buildx-cache - diff --git a/.github/workflows/check.yaml b/.github/workflows/check.yaml new file mode 100644 index 00000000..27bc4d72 --- /dev/null +++ b/.github/workflows/check.yaml @@ -0,0 +1,81 @@ +--- +name: Check πŸ›  + +on: + pull_request: + types: + - opened + - synchronize + - reopened + - ready_for_review + branches: + - main + - develop + push: + branches: + - main + - develop + workflow_dispatch: + +jobs: + audit: + name: Audit Dependencies πŸ•΅οΈβ€β™‚οΈ + uses: insightsengineering/r.pkg.template/.github/workflows/audit.yaml@main + r-cmd: + name: R CMD Check 🧬 + uses: insightsengineering/r.pkg.template/.github/workflows/build-check-install.yaml@main + secrets: + REPO_GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + with: + additional-env-vars: | + CMDSTAN=/root/.cmdstan + CMDSTAN_PATH=/root/.cmdstan + CMDSTANR_NO_VER_CHECK=true + coverage: + name: Coverage πŸ“” + uses: insightsengineering/r.pkg.template/.github/workflows/test-coverage.yaml@main + secrets: + REPO_GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + with: + additional-env-vars: | + CMDSTAN=/root/.cmdstan + CMDSTAN_PATH=/root/.cmdstan + CMDSTANR_NO_VER_CHECK=true + linter: + if: github.event_name == 'pull_request' + name: SuperLinter πŸ¦Έβ€β™€οΈ + uses: insightsengineering/r.pkg.template/.github/workflows/linter.yaml@main + roxygen: + name: Roxygen πŸ…Ύ + uses: insightsengineering/r.pkg.template/.github/workflows/roxygen.yaml@main + secrets: + REPO_GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + with: + auto-update: true + gitleaks: + name: gitleaks πŸ’§ + uses: insightsengineering/r.pkg.template/.github/workflows/gitleaks.yaml@main + spelling: + if: github.event_name == 'pull_request' + name: Spell Check πŸ†Ž + uses: insightsengineering/r.pkg.template/.github/workflows/spelling.yaml@main + links: + if: github.event_name == 'pull_request' + name: Check URLs 🌐 + uses: insightsengineering/r.pkg.template/.github/workflows/links.yaml@main + version: + name: Version Check 🏁 + uses: insightsengineering/r.pkg.template/.github/workflows/version.yaml@main + licenses: + name: License Check πŸƒ + uses: insightsengineering/r.pkg.template/.github/workflows/licenses.yaml@main + style: + if: github.event_name == 'pull_request' + name: Style Check πŸ‘— + uses: insightsengineering/r.pkg.template/.github/workflows/style.yaml@main + with: + auto-update: true + grammar: + if: github.event_name == 'pull_request' + name: Grammar Check πŸ”€ + uses: insightsengineering/r.pkg.template/.github/workflows/grammar.yaml@main diff --git a/.github/workflows/docs.yaml b/.github/workflows/docs.yaml new file mode 100644 index 00000000..ca5313a0 --- /dev/null +++ b/.github/workflows/docs.yaml @@ -0,0 +1,48 @@ +--- +name: Docs πŸ“š + +on: + push: + branches: + - main + - develop + paths: + - "inst/templates/**" + - "_pkgdown.y[a]ml" + - DESCRIPTION + - "**.md" + - "**.Rmd" + - "man/**" + - "LICENSE.*" + - NAMESPACE + pull_request: + types: + - opened + - synchronize + - reopened + - ready_for_review + branches: + - main + - develop + paths: + - "inst/templates/**" + - "_pkgdown.y[a]ml" + - DESCRIPTION + - "**.md" + - "**.Rmd" + - "man/**" + - "LICENSE.*" + - NAMESPACE + workflow_dispatch: + +jobs: + docs: + name: Pkgdown Docs πŸ“š + uses: insightsengineering/r.pkg.template/.github/workflows/pkgdown.yaml@main + secrets: + REPO_GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + with: + additional-env-vars: | + CMDSTAN=/root/.cmdstan + CMDSTAN_PATH=/root/.cmdstan + CMDSTANR_NO_VER_CHECK=true diff --git a/.github/workflows/on_push_pull.yaml b/.github/workflows/on_push_pull.yaml deleted file mode 100755 index 87f877ff..00000000 --- a/.github/workflows/on_push_pull.yaml +++ /dev/null @@ -1,57 +0,0 @@ -on: - push: - branches: [main, master] - pull_request: - branches: [main, master] - -name: on_push_pull - -jobs: - testthat: - runs-on: ubuntu-latest - - name: jmpost-latest-testthat - - container: - image: ghcr.io/genentech/jmpost:latest - - env: - GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} - - steps: - - uses: actions/checkout@v2 - - - name: testthat - run: | - options(crayon.enabled = TRUE, cli.dynamic = FALSE) - devtools::test(stop_on_failure = TRUE, reporter = testthat::CheckReporter) - shell: Rscript {0} - - is-documented: - runs-on: ubuntu-latest - name: check-is-documented - container: - image: ghcr.io/genentech/jmpost:latest - steps: - - - name: Checkout - uses: actions/checkout@v2 - - - name: Document Code - run: | - options(crayon.enabled = TRUE, cli.dynamic = FALSE) - devtools::document() - shell: Rscript {0} - - - name: Check Is Clean - shell: bash - run: | - git config --global --add safe.directory $(pwd) - git status - if [ -z "$(git status --porcelain)" ]; then - echo "Is Clean" - else - echo "Changed Detected" - exit 2 - fi - diff --git a/.github/workflows/release.yaml b/.github/workflows/release.yaml new file mode 100644 index 00000000..2599619a --- /dev/null +++ b/.github/workflows/release.yaml @@ -0,0 +1,39 @@ +--- +name: Release 🎈 + +on: + push: + tags: + - "v*" + workflow_dispatch: + +jobs: + build: + name: Build package 🎁 + needs: release + uses: insightsengineering/r.pkg.template/.github/workflows/build-check-install.yaml@main + secrets: + REPO_GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + with: + skip-r-cmd-check: true + skip-r-cmd-install: true + additional-env-vars: | + CMDSTAN=/root/.cmdstan + CMDSTAN_PATH=/root/.cmdstan + CMDSTANR_NO_VER_CHECK=true + docs: + name: Pkgdown Docs πŸ“š + needs: release + uses: insightsengineering/r.pkg.template/.github/workflows/pkgdown.yaml@main + secrets: + REPO_GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + with: + additional-env-vars: | + CMDSTAN=/root/.cmdstan + CMDSTAN_PATH=/root/.cmdstan + CMDSTANR_NO_VER_CHECK=true + release: + name: Create release πŸŽ‰ + uses: insightsengineering/r.pkg.template/.github/workflows/release.yaml@main + permissions: + contents: write From da7ef2becec38370553f7cb8c55a74888b5e5ddb Mon Sep 17 00:00:00 2001 From: cicdguy <26552821+cicdguy@users.noreply.github.com> Date: Wed, 10 May 2023 06:21:44 -0500 Subject: [PATCH 2/8] Add wordlist --- DESCRIPTION | 3 +- NEWS.md | 3 ++ README.md | 2 +- inst/WORDLIST | 66 ++++++++++++++++++++++++++++++++ inst/stan/base/longitudinal.stan | 2 +- tests/spelling.R | 3 ++ vignettes/lm_gsf.md | 2 +- 7 files changed, 77 insertions(+), 4 deletions(-) create mode 100644 NEWS.md create mode 100644 inst/WORDLIST create mode 100644 tests/spelling.R diff --git a/DESCRIPTION b/DESCRIPTION index 159bfb15..3d40d486 100755 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -9,6 +9,7 @@ Authors@R: c( person("F. Hoffmann-La Roche AG", role = c("cph", "fnd")) ) Description: TODO +Language: en-US License: Apache License (>= 2) Encoding: UTF-8 Roxygen: list(markdown = TRUE) @@ -21,7 +22,7 @@ Imports: statmod, Matrix, jinjar -Suggests: +Suggests: testthat (>= 3.0.0) Config/testthat/edition: 3 Collate: diff --git a/NEWS.md b/NEWS.md new file mode 100644 index 00000000..4985bf94 --- /dev/null +++ b/NEWS.md @@ -0,0 +1,3 @@ +# jmpost 0.0.0.9000 + +* In development diff --git a/README.md b/README.md index 352f9e74..7bf97817 100755 --- a/README.md +++ b/README.md @@ -9,7 +9,7 @@ More specifically, the model implemented in this package utilizes a modeling fra **[1]** [Tardivon _et al._ Association between tumor size kinetics and survival in patients with urothelial carcinoma treated with atezolizumab: Implications for patient follow-up. _Clin Pharm Ther_, 2019](https://doi.org/10.1002/cpt.1450). **[2]** [Kerioui _et al._ Bayesian inference using Hamiltonian Monte-Carlo algorithm for nonlinear joint modeling in the context of cancer immunotherapy. _Stat in Med_, 2020](https://doi.org/10.1002/sim.8756). -**[3]** [Kerioui _et al._ Modelling the association between biomarkers and clinical outcome: An introduction to nonlinear joint models. _Br J Clin Pharm_, 2022](https://doi.org/10.1111/bcp.15200). +**[3]** [Kerioui _et al._ Modeling the association between biomarkers and clinical outcome: An introduction to nonlinear joint models. _Br J Clin Pharm_, 2022](https://doi.org/10.1111/bcp.15200). The models are implemented in [STAN](https://mc-stan.org/), and the package provides a flexible user interface. Please reach out to us via issues or email (see the `DESCRIPTION` file) if you have comments or questions, thank you! diff --git a/inst/WORDLIST b/inst/WORDLIST new file mode 100644 index 00000000..2caa7e7f --- /dev/null +++ b/inst/WORDLIST @@ -0,0 +1,66 @@ +Clin +Kerioui +LogLogistic +Nind +Overal +Pharm +SLD +StanAll +StanMod +StanModule +Tardivon +Ther +al +atezolizumab +biomarker +biomarkers +cX +coeficient +coloured +colours +conection +contails +direcotry +dsld +et +frac +gi +gl +guasian +ie +ig +ij +immunotherapy +initival +inits +lm +ln +longMod +longditudinal +longmodel +mathcal +mesasurements +modlel +odels +oint +os +parametrisation +rajectories +redicting +si +sk +sl +sld +stan +templated +ther +tikzpicture +tikzstyle +timepoint +ttg +tumour +urothelial +urvival +verall +xl +xshift diff --git a/inst/stan/base/longitudinal.stan b/inst/stan/base/longitudinal.stan index 263c418d..bf0617d1 100755 --- a/inst/stan/base/longitudinal.stan +++ b/inst/stan/base/longitudinal.stan @@ -5,7 +5,7 @@ data{ // Source - base/longitudinal.stan // - // Longditudinal data + // Longitudinal data int Nta_total; // Total number of tumor assessments. int Nta_obs_y; // Number of observed tumor assessments (not censored). int Nta_cens_y; // Number of censored tumor assessments (below threshold). diff --git a/tests/spelling.R b/tests/spelling.R new file mode 100644 index 00000000..6713838f --- /dev/null +++ b/tests/spelling.R @@ -0,0 +1,3 @@ +if(requireNamespace('spelling', quietly = TRUE)) + spelling::spell_check_test(vignettes = TRUE, error = FALSE, + skip_on_cran = TRUE) diff --git a/vignettes/lm_gsf.md b/vignettes/lm_gsf.md index 4b53fbb1..5df45b37 100755 --- a/vignettes/lm_gsf.md +++ b/vignettes/lm_gsf.md @@ -1,7 +1,7 @@ # SLD Model Specification -## Longditudinal Model +## Longitudinal Model $$ y_{ij} \sim \mathcal{N}(SLD_{ij}, SLD_{ij}^2 \sigma^2) From 59ef6568f2f45707da0a75c9638cfc360ab48eaa Mon Sep 17 00:00:00 2001 From: github-actions <41898282+github-actions[bot]@users.noreply.github.com> Date: Wed, 10 May 2023 11:23:58 +0000 Subject: [PATCH 3/8] [skip actions] Restyle files --- tests/spelling.R | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/tests/spelling.R b/tests/spelling.R index 6713838f..13f77d96 100644 --- a/tests/spelling.R +++ b/tests/spelling.R @@ -1,3 +1,6 @@ -if(requireNamespace('spelling', quietly = TRUE)) - spelling::spell_check_test(vignettes = TRUE, error = FALSE, - skip_on_cran = TRUE) +if (requireNamespace("spelling", quietly = TRUE)) { + spelling::spell_check_test( + vignettes = TRUE, error = FALSE, + skip_on_cran = TRUE + ) +} From ee5f57b7974ae7d3bb83c52a4ec9ca45b083be72 Mon Sep 17 00:00:00 2001 From: cicdguy <26552821+cicdguy@users.noreply.github.com> Date: Wed, 10 May 2023 08:13:21 -0500 Subject: [PATCH 4/8] Delint markdown files --- README.md | 6 ++++-- tests/spelling.R | 6 ------ vignettes/Intro_to_jmpost.Rmd | 20 ++++++++++---------- vignettes/data_objects.md | 20 +++----------------- vignettes/initial_values.md | 14 ++------------ vignettes/lm_random_slope.md | 3 --- vignettes/models_survival.md | 13 ++++--------- 7 files changed, 23 insertions(+), 59 deletions(-) delete mode 100644 tests/spelling.R diff --git a/README.md b/README.md index 7bf97817..1184c184 100755 --- a/README.md +++ b/README.md @@ -1,14 +1,16 @@ +# jmpost `jmpost` is currently under development and not yet testable. The goal of the `jmpost` package is to fit joint models involving: + 1. a parametric time-to-event sub-model, 2. a nonlinear (or linear) mixed-effect sub-model, describing individual time profiles (_i.e._ trajectories) for a continuous marker, 3. a link function (_a.k.a._ association term). More specifically, the model implemented in this package utilizes a modeling framework described previously **[1-3]** to link overall survival to tumor size data in oncology clinical trials. -**[1]** [Tardivon _et al._ Association between tumor size kinetics and survival in patients with urothelial carcinoma treated with atezolizumab: Implications for patient follow-up. _Clin Pharm Ther_, 2019](https://doi.org/10.1002/cpt.1450). -**[2]** [Kerioui _et al._ Bayesian inference using Hamiltonian Monte-Carlo algorithm for nonlinear joint modeling in the context of cancer immunotherapy. _Stat in Med_, 2020](https://doi.org/10.1002/sim.8756). +**[1]** [Tardivon _et al._ Association between tumor size kinetics and survival in patients with urothelial carcinoma treated with atezolizumab: Implications for patient follow-up. _Clin Pharm Ther_, 2019](https://doi.org/10.1002/cpt.1450). +**[2]** [Kerioui _et al._ Bayesian inference using Hamiltonian Monte-Carlo algorithm for nonlinear joint modeling in the context of cancer immunotherapy. _Stat in Med_, 2020](https://doi.org/10.1002/sim.8756). **[3]** [Kerioui _et al._ Modeling the association between biomarkers and clinical outcome: An introduction to nonlinear joint models. _Br J Clin Pharm_, 2022](https://doi.org/10.1111/bcp.15200). The models are implemented in [STAN](https://mc-stan.org/), and the package provides a flexible user interface. diff --git a/tests/spelling.R b/tests/spelling.R deleted file mode 100644 index 13f77d96..00000000 --- a/tests/spelling.R +++ /dev/null @@ -1,6 +0,0 @@ -if (requireNamespace("spelling", quietly = TRUE)) { - spelling::spell_check_test( - vignettes = TRUE, error = FALSE, - skip_on_cran = TRUE - ) -} diff --git a/vignettes/Intro_to_jmpost.Rmd b/vignettes/Intro_to_jmpost.Rmd index 8e2f29f3..dec36dec 100755 --- a/vignettes/Intro_to_jmpost.Rmd +++ b/vignettes/Intro_to_jmpost.Rmd @@ -6,13 +6,13 @@ author: - Georgios Kazantzidis - Daniel SabanΓ©s BovΓ© - Xiaoyan Fang -header-includes: +header-includes: - \usepackage{tikz} - \usepackage{pgfplots} - \usepackage{graphicx} - \usetikzlibrary{mindmap} - - \usetikzlibrary{positioning} - - \usetikzlibrary{shapes.geometric, arrows} + - \usetikzlibrary{positioning} + - \usetikzlibrary{shapes.geometric, arrows} output: pdf_document: default bibliography: bibliography.bib @@ -20,7 +20,7 @@ vignette: | %\VignetteIndexEntry{Introduction to `jmpost`} %\VignetteEncoding{UTF-8} %\VignetteEngine{knitr::rmarkdown} - + --- ```{r, include = FALSE} @@ -32,14 +32,14 @@ knitr::opts_chunk$set( ``` # Introduction -In clinical development of new cancer medicines a key challenge is to decide on the next steps after having collected the initial data in Phase 1 studies. -In phase 1 studies patients are followed until the cancer is growing again (progression event). -Thus, reliable overall survival data are not available. +In clinical development of new cancer medicines a key challenge is to decide on the next steps after having collected the initial data in Phase 1 studies. +In phase 1 studies patients are followed until the cancer is growing again (progression event). +Thus, reliable overall survival data are not available. However, we can use additional information from previous clinical trials or real-world data to understand the time between progression and death. Specifically, we can correlate the tumor response data, that are longitudinal measurements, (i.e. the shrinkage and growth of tumor lesions over time) with the overall survival of the patients or their hazard ratios [@beyer2020multistate]. Thereby we can predict the overall survival in our Phase 1 studies and therefore make better decisions [@kerioui2020bayesian]. -The sum of the longest diameters (SLD) of the cancer lesions provides an effective and representative biomarker for longitudinal measurements in solid cancer. +The sum of the longest diameters (SLD) of the cancer lesions provides an effective and representative biomarker for longitudinal measurements in solid cancer. SLD is relatively easy to measure with reasonable accuracy while also being a good representative of the development of cancer. The association between the longitudinal data and the overall survival can be achieved in different ways. @@ -51,10 +51,10 @@ Firstly, the use of historical data from previous trials is improving the perfor Secondly, prior definition allows the adjustment of the models according to the specific factors of each analysis. Although these statistical methods provide significant advantages in decision making processes in cancer drug development, their implementation is relatively intricate and time consuming. -Compared to other R-packages such as `joineR` and `jmbayes`, `jmpost` provides the flexibility of non-linear mixed effects modeling for the longitudinal part of the joint model. +Compared to other R-packages such as `joineR` and `jmbayes`, `jmpost` provides the flexibility of non-linear mixed effects modeling for the longitudinal part of the joint model. Non-linear parametrization respects the biological implications of clinical trials in oncology where the observed tumor response usually follows non-linear patterns. -# Minimal usage +# Minimal usage The minimal usage of the `jmpost` package calls example stan code already provided in the repository of the package. The following chunk illustrates the absolute minimal usage of the package. diff --git a/vignettes/data_objects.md b/vignettes/data_objects.md index 0d0ba719..4e4642d1 100644 --- a/vignettes/data_objects.md +++ b/vignettes/data_objects.md @@ -1,7 +1,4 @@ - - -## Global Data Objects - +# Global Data Objects ```R int Nind # Number of individuals @@ -24,11 +21,8 @@ array[n_arms] int n_index_per_arm array[Nind] int index_per_arm ``` - - ## Survival Data Objects - ```R int Nind_dead # Number of patients who died vector[Nind] Times # Event and Censor times per subject @@ -42,10 +36,8 @@ vector[n_nodes] nodes vector[n_nodes] weights ``` - ## Longitudinal Data Objects - ```R int Nta_total # Number of tumour observations int Nta_obs_y # Number of observed tumour values @@ -53,14 +45,14 @@ int Nta_cens_y # Number of censored tumour values # Nta_total = Nta_obs_y + Nta_cens_y -# Patient index numbers for ownership of tumour observations +# Patient index numbers for ownership of tumour observations array[Nta_total] int ind_index # Tumour Index numbers for observed tumour values array[Nta_obs_y] int obs_y_index # Tumour Index numbers for censored tumour values -array[Nta_cens_y] int cens_y_index +array[Nta_cens_y] int cens_y_index vector[Nta_total] Yobs # Tumour size observations @@ -86,10 +78,4 @@ array [3] int n_mat_inds_cens_y; vector[n_mat_inds_cens_y[1]] w_mat_inds_cens_y; array[n_mat_inds_cens_y[2]] int v_mat_inds_cens_y; array[n_mat_inds_cens_y[3]] int u_mat_inds_cens_y; - - - ``` - - - diff --git a/vignettes/initial_values.md b/vignettes/initial_values.md index e08e91e8..74f41678 100644 --- a/vignettes/initial_values.md +++ b/vignettes/initial_values.md @@ -3,9 +3,9 @@ You can specify initial values in one of three ways: -1. Use the default. By default all initial values are set to the mean of the prior distribution. That is if you use a `prior_gamma(1, 3)` the initial value will be set to `1/3`. +1. Use the default. By default all initial values are set to the mean of the prior distribution. That is if you use a `prior_gamma(1, 3)` the initial value will be set to `1/3`. -2. Manually specify in the prior function. That is you can set initial values via `prior_normal(0, 1, init = 20)`. Note that if the parameter represents multiple parameters (e.g. 1 per arm) then single values are replicated for each of the parameters. You can specify different initival values by instead specifying a vector e.g. `prior_normal(0, 1, init = c(10, 20))`. Care is needed to ensure the correct dimensionality. +2. Manually specify in the prior function. That is you can set initial values via `prior_normal(0, 1, init = 20)`. Note that if the parameter represents multiple parameters (e.g. 1 per arm) then single values are replicated for each of the parameters. You can specify different initival values by instead specifying a vector e.g. `prior_normal(0, 1, init = c(10, 20))`. Care is needed to ensure the correct dimensionality. 3. Provide a list of initial values to the sampling function. Instead of relying on `jmpost` to generate the initial values for you, you can instead specify a list of values directly to the sampling function e.g. @@ -30,13 +30,3 @@ samples <- mp <- sampleStanModel( init = initial_values ) ``` - - - - - - - - - - diff --git a/vignettes/lm_random_slope.md b/vignettes/lm_random_slope.md index 89901425..fbae6e13 100755 --- a/vignettes/lm_random_slope.md +++ b/vignettes/lm_random_slope.md @@ -1,8 +1,6 @@ # Longitudinal Model Specifications - - ## Random Slope Model $$ @@ -31,4 +29,3 @@ Where: - $\phi$ is a global scaling coeficient That is to say a subjects hazard is constant where that constant is a multiple of their random slope - diff --git a/vignettes/models_survival.md b/vignettes/models_survival.md index a526db8b..0a0da69c 100755 --- a/vignettes/models_survival.md +++ b/vignettes/models_survival.md @@ -1,12 +1,9 @@ - - - # Overal Survival Model Specification ## General Overal Survival Model $$ -h_i(t) = +h_i(t) = h_{0}(t | \theta) exp\left( \beta X_i + C(t \mid \phi_i) @@ -19,8 +16,6 @@ where: * $X_i$ is the matrix of covariates for subject $i$ * $C(t \mid \phi_i)$ is the link contribution function at time $t$ given parameters $\phi_i$ from a corresponding longditudinal - - ## Weibull Distribution @@ -33,15 +28,15 @@ $$ $$ \begin{aligned} -t &\sim \text{Log-Logistic}(\lambda, p) +t &\sim \text{Log-Logistic}(\lambda, p) \\ \\ f_0(t \mid \lambda, p) &= \frac{ \lambda p (\lambda t)^{p-1} } { (1 + (\lambda t)^p)^2 -} -\\ +} +\\ \\ h_0(t \mid \lambda, p) &= \frac{ \lambda p (\lambda t) ^{p -1} From 95e27b0b3eacef99ba7b68e19f5c4706ab2df598 Mon Sep 17 00:00:00 2001 From: github-actions <41898282+github-actions[bot]@users.noreply.github.com> Date: Wed, 10 May 2023 13:23:06 +0000 Subject: [PATCH 5/8] [skip actions] Restyle files --- vignettes/Intro_to_jmpost.Rmd | 68 ++++++++++++++++++----------------- 1 file changed, 35 insertions(+), 33 deletions(-) diff --git a/vignettes/Intro_to_jmpost.Rmd b/vignettes/Intro_to_jmpost.Rmd index dec36dec..e3f03297 100755 --- a/vignettes/Intro_to_jmpost.Rmd +++ b/vignettes/Intro_to_jmpost.Rmd @@ -66,7 +66,7 @@ my_long_mod <- LongMod() my_temp_os <- LogLogisticOs() # The link between longitudinal and OS models my_link <- ttgLink() -# The complete survival submodel +# The complete survival submodel my_os <- parametrize(osmod = my_temp_os, link = my_link) # The full model JointModel(long = my_long_mod, os = my_os) @@ -120,16 +120,16 @@ The longitudinal sub-model, including all the required information for the model ```{r Longitudinal model} Long_mod <- LongModel( - stan = StanModule( - functions = "exp_long_functions.stan", - data = "exp_long_data.stan", - priors = long_prior(), - model = "", - inits = list(), - parameters = "exp_long_parameters.stan", - transformed_parameters = "exp_long_transformed_parametes.stan", - generated_quantities = "" - ) + stan = StanModule( + functions = "exp_long_functions.stan", + data = "exp_long_data.stan", + priors = long_prior(), + model = "", + inits = list(), + parameters = "exp_long_parameters.stan", + transformed_parameters = "exp_long_transformed_parametes.stan", + generated_quantities = "" + ) ) ``` @@ -141,10 +141,10 @@ The user is able to populate the different parts of the object with stan code ei ```{r Longitudinal model2 } example_stan <- LongModel( - stan = StanModule( - functions = "type stan code here;", - data = "type stan code here;" - ) + stan = StanModule( + functions = "type stan code here;", + data = "type stan code here;" + ) ) ``` @@ -187,14 +187,15 @@ Figure~2 depicts the structure of a `StanModule` object. Gray coloured cells indicate character strings containing stan code whereas yellow coloured cells mark list objects. ```{r} -StanModule(functions = "", - data = "", - model = "", - parameters = "", - transformed_parameters = "", - generated_quantities = "", - priors = list(), - inits = list() +StanModule( + functions = "", + data = "", + model = "", + parameters = "", + transformed_parameters = "", + generated_quantities = "", + priors = list(), + inits = list() ) ``` @@ -242,9 +243,9 @@ The link is an object containing the parametrisation of the longitudinal sub-mod ```{r} link <- HazardLink( - parameters = "beta_ttg", - contribution = "* rep_matrix(ttg(psi_ks, psi_kg, psi_phi), rows(time))", - stan = StanModule( ) + parameters = "beta_ttg", + contribution = "* rep_matrix(ttg(psi_ks, psi_kg, psi_phi), rows(time))", + stan = StanModule() ) ``` @@ -273,13 +274,14 @@ jm <- JointModel(long = Long_mod, os = os) Although the stan code of the model is obtained in `jm`, additional information for the model run should be provided. ```{r} -jm_options <- mcmc_options(chains = 1, - parallel_chains = 1, - iter_warmup = 200, - iter_sampling = 500, - max_treedepth = 12, - adapt_delta = .9, - gauss_legendre = gauss_legendre() +jm_options <- mcmc_options( + chains = 1, + parallel_chains = 1, + iter_warmup = 200, + iter_sampling = 500, + max_treedepth = 12, + adapt_delta = .9, + gauss_legendre = gauss_legendre() ) ``` From 5bebc32a9bb11c43e7017af97725be95429b003b Mon Sep 17 00:00:00 2001 From: cicdguy <26552821+cicdguy@users.noreply.github.com> Date: Wed, 10 May 2023 08:32:05 -0500 Subject: [PATCH 6/8] More delinting --- vignettes/Intro_to_jmpost.Rmd | 123 ++++++++++++++++------------------ vignettes/lm_gsf.md | 30 ++------- vignettes/models_survival.md | 2 - 3 files changed, 65 insertions(+), 90 deletions(-) diff --git a/vignettes/Intro_to_jmpost.Rmd b/vignettes/Intro_to_jmpost.Rmd index e3f03297..3f0bb327 100755 --- a/vignettes/Intro_to_jmpost.Rmd +++ b/vignettes/Intro_to_jmpost.Rmd @@ -27,20 +27,21 @@ vignette: | knitr::opts_chunk$set( collapse = TRUE, comment = "#>", - eval = F + eval = FALSE ) ``` + # Introduction In clinical development of new cancer medicines a key challenge is to decide on the next steps after having collected the initial data in Phase 1 studies. In phase 1 studies patients are followed until the cancer is growing again (progression event). Thus, reliable overall survival data are not available. However, we can use additional information from previous clinical trials or real-world data to understand the time between progression and death. -Specifically, we can correlate the tumor response data, that are longitudinal measurements, (i.e. the shrinkage and growth of tumor lesions over time) with the overall survival of the patients or their hazard ratios [@beyer2020multistate]. +Specifically, we can correlate the tumor response data, that are longitudinal measurements, (i.e. the shrinkage and growth of tumor lesions over time) with the overall survival of the patients or their hazard ratios [@beyer2020multistate]. Thereby we can predict the overall survival in our Phase 1 studies and therefore make better decisions [@kerioui2020bayesian]. The sum of the longest diameters (SLD) of the cancer lesions provides an effective and representative biomarker for longitudinal measurements in solid cancer. -SLD is relatively easy to measure with reasonable accuracy while also being a good representative of the development of cancer. +SLD is relatively easy to measure with reasonable accuracy while also being a good representative of the development of cancer. The association between the longitudinal data and the overall survival can be achieved in different ways. The predicted SLD value or its first derivative over time can be directly treated as the linking factor in the hazard model for survival. @@ -56,7 +57,7 @@ Non-linear parametrization respects the biological implications of clinical tria # Minimal usage -The minimal usage of the `jmpost` package calls example stan code already provided in the repository of the package. +The minimal usage of the `jmpost` package calls example stan code already provided in the repository of the package. The following chunk illustrates the absolute minimal usage of the package. ```{r, eval=FALSE} @@ -73,14 +74,14 @@ JointModel(long = my_long_mod, os = my_os) ``` The final two steps of the package workflow are not available yet: -First, the running step of the model, where the user calls the defined model along with the data and the preferred options. +First, the running step of the model, where the user calls the defined model along with the data and the preferred options. Secondly, the summary report of the model run. # Customizing the analysis -`jmpost` is a package aiming to optimize the procedure of developing a new Bayesian joint model. -The package provides complete examples of joint modeling code. -Here we explain the workflow of the package `jmpost` using `stan` libraries as an example. +`jmpost` is a package aiming to optimize the procedure of developing a new Bayesian joint model. +The package provides complete examples of joint modeling code. +Here we explain the workflow of the package `jmpost` using `stan` libraries as an example. \begin{figure}[!] \tikzstyle{Long} = [rectangle, rounded corners, minimum width=3cm, minimum height=1cm, text centered, draw=black, fill=green!30 ] @@ -93,7 +94,6 @@ Here we explain the workflow of the package `jmpost` using `stan` libraries as a \tikzstyle{arrow} = [thick,->,>=stealth ] - \begin{tikzpicture}[node distance=2cm ] \node (long) [Long ] {Longitudinal}; \node (os_temp) [Os, right of = long, xshift = 4cm ] {Overall survival (templated)}; @@ -101,7 +101,6 @@ Here we explain the workflow of the package `jmpost` using `stan` libraries as a \node (os) [Os, below = 2cm of os_temp] {Overall survival}; \node (model) [jmpost, below left = 2cm and 0.5cm of os] {Full joint model}; - \draw [arrow] (long) -- (link); \draw [arrow] (os_temp) -- (os); \draw [arrow] (link) -- (os); @@ -113,7 +112,6 @@ Here we explain the workflow of the package `jmpost` using `stan` libraries as a \caption{Package workflow. Different colours signify objects of different classes. Arrows indicate the relationships between the objects.} \end{figure} - Figure~1 illustrates the main workflow of the `jmpost` package. Three objects need to be combined for the construction of the full Bayesian joint model. The longitudinal sub-model, including all the required information for the modeling of the longitudinal measurements, the survival sub-model, including information of for the time to event data and the link or association factor, containing information for the conection of the two sub-models. @@ -134,10 +132,10 @@ Long_mod <- LongModel( ``` The construction of the different parts of the model is flexible. -Here, we create the longitudinal part by calling the `LongModel` constructor. -`LongModel` identifies the type of the sub-model whereas the helper constructor `StanModule` is responsible for the actual stan code. +Here, we create the longitudinal part by calling the `LongModel` constructor. +`LongModel` identifies the type of the sub-model whereas the helper constructor `StanModule` is responsible for the actual stan code. The created object contains six compartments containing different parts of the stan code (functions, data, priors, model, parameters, transformed parameters and generated quantities) while also `inits` which contains the initial values for the Markov Chain Monte Carlo run. -The user is able to populate the different parts of the object with stan code either by providing the names of the stan files contained in the `jmpost` directory, the full paths of stan files outside the `jmpost` direcotry or by providing stan code in text form (code chunk below). +The user is able to populate the different parts of the object with stan code either by providing the names of the stan files contained in the `jmpost` directory, the full paths of stan files outside the `jmpost` direcotry or by providing stan code in text form (code chunk below). ```{r Longitudinal model2 } example_stan <- LongModel( @@ -148,7 +146,6 @@ example_stan <- LongModel( ) ``` - \begin{figure} \tikzstyle{Char} = [rectangle, rounded corners, minimum width=3cm, minimum height=1cm, text centered, draw=black, fill=gray!30 ] @@ -183,8 +180,8 @@ example_stan <- LongModel( \caption{Structure of `StanModule` object.} \end{figure} -Figure~2 depicts the structure of a `StanModule` object. -Gray coloured cells indicate character strings containing stan code whereas yellow coloured cells mark list objects. +Figure~2 depicts the structure of a `StanModule` object. +Gray coloured cells indicate character strings containing stan code whereas yellow coloured cells mark list objects. ```{r} StanModule( @@ -200,46 +197,42 @@ StanModule( ``` -The `jmpost` repository contains stan code for the creation of a specific stan model. -The user is able to populate the longitudinal object `Long_mod` either by specifying the stan files containing the relevant code of the section or by directly typing the stan code into R. -We split the model section of the stan file into two parts, the `priors` and the `model`. +The `jmpost` repository contains stan code for the creation of a specific stan model. +The user is able to populate the longitudinal object `Long_mod` either by specifying the stan files containing the relevant code of the section or by directly typing the stan code into R. +We split the model section of the stan file into two parts, the `priors` and the `model`. `Priors` contain a list with the names of the parameters and their prior distribution definitions while `model` is a character that contains model relevant information such us the definition of the likelihood. -The list of priors should be a named list of strings where each element represents the stan code for the definition of the desired density. +The list of priors should be a named list of strings where each element represents the stan code for the definition of the desired density. ```{r} prior_set <- StanModule(priors = list("param_A" = "normal(0,1);")) ``` -The user can modify the prior list at any point during the workflow. +The user can modify the prior list at any point during the workflow. ```{r} prior_set@priors$param_A <- "normal(1,1);" ``` - Call of the object prints the stan code of all the sections of the object. ```{r Long_model print} Long_mod ``` -The overall survival object contails all ther information for the parametrisation of the second sub-model of the joint model. +The overall survival object contails all ther information for the parametrisation of the second sub-model of the joint model. The overall survival object has identical structure as the Longitudinal object (Figure~2). -The slots `functions`, `parameters`, `model`, `transformed parameters` and `generated quantities` should be provided with stan code whereas `priors` and `initial values` should be filled with lists. -An example with Log-logistic parametrisation of the overall survival sub-modlel object is provided in `LogLogisticOS()`. - - +The slots `functions`, `parameters`, `model`, `transformed parameters` and `generated quantities` should be provided with stan code whereas `priors` and `initial values` should be filled with lists. +An example with Log-logistic parametrisation of the overall survival sub-modlel object is provided in `LogLogisticOS()`. ```{r} temp_os <- LogLogisticOs() ``` -Although the Longitudinal and the Survival object identical structures, the contained stan code has some fundamental differences. +Although the Longitudinal and the Survival object identical structures, the contained stan code has some fundamental differences. The functions contained in the overall survival object depend on the parametrisation of the longitudinal part. -Thus, in order to complete the stan code of the former object, additional information is required. -The required details are provided to the overall survival object through the link object. -The link is an object containing the parametrisation of the longitudinal sub-model, that will be used as part of the overall survival object. - +Thus, in order to complete the stan code of the former object, additional information is required. +The required details are provided to the overall survival object through the link object. +The link is an object containing the parametrisation of the longitudinal sub-model, that will be used as part of the overall survival object. ```{r} link <- HazardLink( @@ -249,29 +242,29 @@ link <- HazardLink( ) ``` -Link objects or association factor, contain the slots of `parameters`, `contribution` and `stan`. -Here, the slot `parameters` is the name of the linking parameter of the two sub-models. -`Contribution` contains the stan code relevant to the linking parameter. -`Parameters` and `contribution` contain all the information required for the completion of the overall survival model. +Link objects or association factor, contain the slots of `parameters`, `contribution` and `stan`. +Here, the slot `parameters` is the name of the linking parameter of the two sub-models. +`Contribution` contains the stan code relevant to the linking parameter. +`Parameters` and `contribution` contain all the information required for the completion of the overall survival model. -The merge of the `Longitudinal`, `Overall survival` and `link` object takes place in two parts. -First, the parameterisation of the overall survival sub-model needs to be completed using the information from the link object (part 1). +The merge of the `Longitudinal`, `Overall survival` and `link` object takes place in two parts. +First, the parameterisation of the overall survival sub-model needs to be completed using the information from the link object (part 1). ```{r} os <- parametrize(osmod = temp_os, link = link) ``` -Secondly, the `Longitudinal` and `Survival` sub-models need to be combined (part 2). -Function `JointModel` merges the two sub-models, creates a character containing the full stan model, creates a stan file with the produced model and checks the functionality of the produced file. +Secondly, the `Longitudinal` and `Survival` sub-models need to be combined (part 2). +Function `JointModel` merges the two sub-models, creates a character containing the full stan model, creates a stan file with the produced model and checks the functionality of the produced file. Mistakes in the stan file will be reported from the `JointModel` call. ```{r} jm <- JointModel(long = Long_mod, os = os) ``` -`jm` object contains all the complete stan code for the joint model. +`jm` object contains all the complete stan code for the joint model. `jmpost` uses the package `cmdstanr` for the compilation of the stan model. -Although the stan code of the model is obtained in `jm`, additional information for the model run should be provided. +Although the stan code of the model is obtained in `jm`, additional information for the model run should be provided. ```{r} jm_options <- mcmc_options( @@ -285,23 +278,23 @@ jm_options <- mcmc_options( ) ``` -`mcmc_options` object contains technical information relevant to the model run. -Most of the options are borrowed from the package `cmdstanr` except from the `gauss_legendre` that defines the integration parameters. +`mcmc_options` object contains technical information relevant to the model run. +Most of the options are borrowed from the package `cmdstanr` except from the `gauss_legendre` that defines the integration parameters. -## Example +## Example -### SLD model +### SLD model \begin{align} y_{ij} \sim N(SLD_{ij}, SLD^2_{ij}\sigma^2) \end{align} -Where: +Where: \begin{itemize} \item $y_{ij}$ is the observed tumor mesasurements \item $SLD_{ij}$ is the expected sum of longest diameter for subject $i$ at time point $j$ \end{itemize} -### Expected SLD +### Expected SLD --- @@ -310,14 +303,14 @@ SLD_{ij} = b_i[\phi_ie^{-s_it_{ij}}+(1-\phi_i)e^{g_it_{ij}}] \end{align} Where: \begin{itemize} -\item $i$ is the subject index -\item $j$ is the visit index -\item $t_{ij}$ the time since first treatment for subject $i$ at visit $j$ +\item $i$ is the subject index +\item $j$ is the visit index +\item $t_{ij}$ the time since first treatment for subject $i$ at visit $j$ \item $SLD_{ij}$ is the observed SLD measurement for subject $i$ at visit $j$ \item $b_i$ is the Baseline sld measurement -\item $s_i$ is the kinetics shrinkage parameter -\item $g_i$ is the kinetics tumor growth parameter -\item $\phi_i$ is the proportion of cells affected by the treatment +\item $s_i$ is the kinetics shrinkage parameter +\item $g_i$ is the kinetics tumor growth parameter +\item $\phi_i$ is the proportion of cells affected by the treatment \end{itemize} ### The SLD parameters @@ -333,11 +326,11 @@ b_i = exp(ln(m_{bl_i}) + \eta_{bi}*\omega_b) Where: \begin{itemize} -\item $i$ is the subject index +\item $i$ is the subject index \item $l_i$ is the group/treatment index for subject $i$ \item $m_{xl_i}$ is the mean for parameter $x$ in group $l_i$ \item $\eta_{xi}$ is a random effects offset on parameter $x$ for subject $i$ -\item $\omega_x$ is the variance for the random effects on parameter $x$ +\item $\omega_x$ is the variance for the random effects on parameter $x$ \end{itemize} ### The SLD Hyper-parameters @@ -353,11 +346,11 @@ logit(m_{\phi} l_i) \sim N(logit(\mu_{\phi}), \sigma_{\phi}) Where: \begin{itemize} \item $l_i$ is the group/treatment index of subject $i$ -\item $\mu_x$ is the mean of the parameter distribution +\item $\mu_x$ is the mean of the parameter distribution \item $\sigma_x$ is the variance of the parameter distribution \end{itemize} -### Survival model +### Survival model --- @@ -372,7 +365,7 @@ Where: \item $h_0(t)$ is the baseline hazard function \item $f_i(t)$ is the derivative of the SLD trajectory for subject $i$ at time $t$ \item $G_i$ is the time to growth for subject $i$ -\item $X_i$ is the matrix of covariates at each visit for subject $i$ +\item $X_i$ is the matrix of covariates at each visit for subject $i$ \item $\beta, \gamma \& \beta_c$ are the strength coefficients for the derivative of the SLD trajectory, the time to growth and the subjects covariates respectively \end{itemize} @@ -389,12 +382,12 @@ S_0(t) = \frac{1}{1+(\lambda t)^p} Where: \begin{itemize} -\item $f(t)$ is the probability density function -\item $h_0(t)$ is the hazard distribution function +\item $f(t)$ is the probability density function +\item $h_0(t)$ is the hazard distribution function \item $S_0(t)$ is the survival distribution function \end{itemize} -### Time to growth +### Time to growth --- @@ -405,8 +398,8 @@ G_i = \frac{logit(\phi_i)+log(s_ig_i^{-1})}{s_i+g_i} Where: \begin{itemize} \item $G_i$ is the time to growth for subject $i$ -\item $b_i$ is the baseline SLD measurement -\item $s_i$ is the kinetics shrinkage parameter +\item $b_i$ is the baseline SLD measurement +\item $s_i$ is the kinetics shrinkage parameter \item $g_i$ is is the kinetics tumor growth parameter \item $\phi_i$ is the proportion of cells affected by the treatment \end{itemize} diff --git a/vignettes/lm_gsf.md b/vignettes/lm_gsf.md index 5df45b37..0109968a 100755 --- a/vignettes/lm_gsf.md +++ b/vignettes/lm_gsf.md @@ -1,6 +1,5 @@ # SLD Model Specification - ## Longitudinal Model $$ @@ -12,7 +11,6 @@ Where: * $y_{ij}$ is the observed tumour mesasurements * $SLD_{ij}$ is the expected sum of longest diameter for subject $i$ at time point $j$ - ### Expected SLD $$ @@ -23,7 +21,7 @@ SLD_{ij} =b_{i} \right] $$ -Where: +Where: * $i$ is the subject index * $j$ is the visit index @@ -34,22 +32,21 @@ Where: * $g_i$ is the kinetics tumour growth parameter * $\phi_i$ is the proportion of cells affected by the treatment - ### Expected SLD Parameters $$ \begin{aligned} \phi_i &= \text{logit}^{-1}\left( \text{logit}(m_{\phi l_i}) + \eta_{\phi i} * \omega_{\phi} -\right) +\right) \\ s_i &= \exp\left( \ln(m_{s l_i}) + \eta_{s i} * \omega_{s} -\right) +\right) \\ g_i &= \exp\left( \ln(m_{g l_i}) + \eta_{g i} * \omega_{g} -\right) +\right) \\ b_i &= \exp\left( \ln(m_{b l_i}) + \eta_{b i} * \omega_{b} @@ -57,7 +54,6 @@ b_i &= \exp\left( \end{aligned} $$ - Where: * $i$ is the subject index @@ -66,7 +62,6 @@ Where: * $\eta_{xi}$ is a random effects offset on parameter $x$ for subject $i$ * $\omega_{x}$ is the variance for the random effects on parameter $x$ - ### Expected SLD Hyper-parameters @@ -84,7 +79,6 @@ Where: * $\mu_x$ is the mean of the parameter distribution * $\sigma_x$ is the variance of the parameter distribution - ### Priors - TODO - remove as not fixed $$ @@ -102,7 +96,7 @@ $$ \mu_{bi} &\sim \text{Lognormal}(55,5); \\ \\ \mu_s &\sim \text{Lognormal}(1,0.5); \\ -\mu_g &\sim \text{Lognormal}(-0.36,1); \\ +\mu_g &\sim \text{Lognormal}(-0.36,1); \\ \mu_{\phi} &\sim Beta(5,5); \\ \\ \sigma_s &\sim \text{Lognormal}(0,0.5); \\ @@ -111,7 +105,6 @@ $$ \end{aligned} $$ - ## Link Contribution $$ @@ -127,16 +120,12 @@ Where: * $G_k(.)$ are arbitrary functions of the SLD parameters. Several available functions provided by this package are listed below * $\tau_k$ is a global scaling coeficient for $G_k(.)$ - - - ### Derivative of the SLD Trajectory (dsld-link) - $$ -G(t \mid b_i, s_i, g_i, \phi_i) = b_i +G(t \mid b_i, s_i, g_i, \phi_i) = b_i \left[ - (1-\phi_i) g_i e^{g_i t} - + (1-\phi_i) g_i e^{g_i t} - \phi_i s_i e^{-s_i t} \right] $$ @@ -150,10 +139,8 @@ Where: (See the "SLD Model" section for full details) - ### Time to Growth (ttg-link) - $$ G(t \mid b_i, s_i, g_i, \phi_i) = \frac{ \text{logit}(\phi_i) + log(s_i g_i^{-1}) @@ -169,6 +156,3 @@ Where: * $s_i$ is the kinetics shrinkage parameter * $g_i$ is the kinetics tumour growth parameter * $\phi_i$ is the proportion of cells affected by the treatment - - - diff --git a/vignettes/models_survival.md b/vignettes/models_survival.md index 0a0da69c..9e09ef01 100755 --- a/vignettes/models_survival.md +++ b/vignettes/models_survival.md @@ -18,12 +18,10 @@ where: ## Weibull Distribution - $$ h_{0}(t \mid \lambda, \gamma) = \lambda \gamma t^{(\gamma - 1)} $$ - ## Log-Logistic Distribution $$ From abe1da5ee2ece3e92dde0ec280ca4dba420b6145 Mon Sep 17 00:00:00 2001 From: cicdguy <26552821+cicdguy@users.noreply.github.com> Date: Wed, 10 May 2023 08:33:52 -0500 Subject: [PATCH 7/8] Restyle package --- R/DataJoint.R | 65 ++-- R/DataLongitudinal.R | 253 ++++++++------- R/DataSurvival.R | 307 +++++++++--------- R/JointModel.R | 185 +++++------ R/Link.R | 52 ++- R/LinkNone.R | 20 +- R/LongitudinalGSF.R | 204 ++++++------ R/LongitudinalModel.R | 31 +- R/LongitudinalRandomSlope.R | 77 ++--- R/Parameter.R | 88 +++--- R/ParameterList.R | 145 +++++---- R/Prior.R | 110 ++++--- R/StanModel.R | 43 ++- R/StanModule.R | 352 ++++++++++----------- R/SurvivalExponential.R | 26 +- R/SurvivalLoglogistic.R | 37 +-- R/SurvivalModel.R | 27 +- R/SurvivalWeibullPH.R | 27 +- R/defaults.R | 72 ++--- R/generics.R | 47 ++- R/simulations.R | 422 ++++++++++++------------- R/utilities.R | 178 +++++------ misc/example_of_use.R | 132 ++++---- misc/tests/gsf-no-link.R | 198 ++++++------ misc/tests/gsf-with-link.R | 164 +++++----- misc/tests/os-exponential.R | 84 +++-- misc/tests/os-loglogistic.R | 81 +++-- misc/tests/os-weibull.R | 85 +++-- misc/tests/rs-no-link.R | 126 ++++---- misc/tests/rs-with-link.R | 128 ++++---- misc/tests/stan_functions.R | 123 ++++--- tests/testthat/test-DataJoint.R | 170 +++++----- tests/testthat/test-DataLongitudinal.R | 54 ++-- tests/testthat/test-DataSurvival.R | 94 +++--- tests/testthat/test-JointModel.R | 28 +- tests/testthat/test-LinkNone.R | 29 +- tests/testthat/test-LinkRandomSlope.R | 42 +-- tests/testthat/test-ParameterList.R | 46 ++- tests/testthat/test-Parameters.R | 27 +- tests/testthat/test-Prior.R | 20 +- tests/testthat/test-compile.R | 16 +- tests/testthat/test-utilities.R | 118 ++++--- 42 files changed, 2145 insertions(+), 2388 deletions(-) diff --git a/R/DataJoint.R b/R/DataJoint.R index 862b9d02..26855fbf 100755 --- a/R/DataJoint.R +++ b/R/DataJoint.R @@ -4,11 +4,11 @@ NULL #' @rdname DataJoint .DataJoint <- setClass( - Class = "DataJoint", - representation = list( - survival = "DataSurvival", - longitudinal = "DataLongitudinal" - ) + Class = "DataJoint", + representation = list( + survival = "DataSurvival", + longitudinal = "DataLongitudinal" + ) ) @@ -31,48 +31,48 @@ NULL #' @export #' @aliases as.list,DataJoint-method DataJoint <- function(survival, longitudinal) { - .DataJoint( - survival = survival, - longitudinal = longitudinal - ) + .DataJoint( + survival = survival, + longitudinal = longitudinal + ) } setValidity( - Class = "DataJoint", - method = function(object) { - lm <- as.data.frame(object@longitudinal) - lvars <- extractVariableNames(object@longitudinal) - os <- as.data.frame(object@survival) - ovars <- extractVariableNames(object@survival) + Class = "DataJoint", + method = function(object) { + lm <- as.data.frame(object@longitudinal) + lvars <- extractVariableNames(object@longitudinal) + os <- as.data.frame(object@survival) + ovars <- extractVariableNames(object@survival) - if (!all(as.character(lm[[lvars$pt]]) %in% as.character(os[[ovars$pt]]))) { - return("There are subjects in the longitudinal data that do not exist in the survival data") - } - if (!all(as.character(os[[ovars$pt]]) %in% as.character(lm[[lvars$pt]]))) { - return("There are subjects in the survival data that do not exist in the longitudinal data") - } + if (!all(as.character(lm[[lvars$pt]]) %in% as.character(os[[ovars$pt]]))) { + return("There are subjects in the longitudinal data that do not exist in the survival data") + } + if (!all(as.character(os[[ovars$pt]]) %in% as.character(lm[[lvars$pt]]))) { + return("There are subjects in the survival data that do not exist in the longitudinal data") } + } ) #' @export setMethod( - "as.list", - signature = "DataJoint", - definition = function(x) { - append( - as.list(x@survival), - as.list(x@longitudinal) - ) - } + "as.list", + signature = "DataJoint", + definition = function(x) { + append( + as.list(x@survival), + as.list(x@longitudinal) + ) + } ) #' @export setAs( - from = "DataJoint", - to = "list", - def = function(from) as.list(from) + from = "DataJoint", + to = "list", + def = function(from) as.list(from) ) @@ -99,4 +99,3 @@ setAs( # which(osd_final$ARM == "MOR_ASG"), # which(osd_final$ARM == "MOR_AEV") # ), - diff --git a/R/DataLongitudinal.R b/R/DataLongitudinal.R index 1bd077f5..ae13cec0 100644 --- a/R/DataLongitudinal.R +++ b/R/DataLongitudinal.R @@ -1,4 +1,3 @@ - #' @include generics.R NULL @@ -7,13 +6,13 @@ setClassUnion("numeric_or_NULL", c("numeric", "NULL")) #' @rdname DataLongitudinal .DataLongitudinal <- setClass( - Class = "DataLongitudinal", - representation = list( - data = "data.frame", - formula = "formula", - subject = "character", - threshold = "numeric_or_NULL" - ) + Class = "DataLongitudinal", + representation = list( + data = "data.frame", + formula = "formula", + subject = "character", + threshold = "numeric_or_NULL" + ) ) #' `DataLongitudinal` @@ -37,147 +36,143 @@ setClassUnion("numeric_or_NULL", c("numeric", "NULL")) #' #' @export DataLongitudinal <- function(data, formula, subject, threshold = NULL) { - .DataLongitudinal( - data = remove_missing_rows(data, formula, c(subject)), - formula = formula, - subject = subject, - threshold = threshold - ) + .DataLongitudinal( + data = remove_missing_rows(data, formula, c(subject)), + formula = formula, + subject = subject, + threshold = threshold + ) } setValidity( - "DataLongitudinal", - function(object) { - if (length(object@subject) > 1) { - return("`subject` should be a length 1 character vector or NULL") - } - if (!length(object@formula) == 3) { - return("`formula` should be a 2 sided formula") - } - if (!length(object@formula[[3]]) == 1) { - return("the RHS of `formula` should only have 1 value") - } - if (! object@subject %in% names(object@data)) { - return("`subject` does not exist in `data`") - } - pt <- object@data[[object@subject]] - if (!(is(pt, "character") | is(pt, "factor"))) { - return("`data[[subject]]` should be of type character or factor") - } - return(TRUE) + "DataLongitudinal", + function(object) { + if (length(object@subject) > 1) { + return("`subject` should be a length 1 character vector or NULL") + } + if (!length(object@formula) == 3) { + return("`formula` should be a 2 sided formula") + } + if (!length(object@formula[[3]]) == 1) { + return("the RHS of `formula` should only have 1 value") + } + if (!object@subject %in% names(object@data)) { + return("`subject` does not exist in `data`") + } + pt <- object@data[[object@subject]] + if (!(is(pt, "character") | is(pt, "factor"))) { + return("`data[[subject]]` should be of type character or factor") } + return(TRUE) + } ) setMethod( - "as.data.frame", - signature = "DataLongitudinal", - definition = function(x, row.names, optional, ...) { - df_fct <- x@data - vars <- extractVariableNames(x) - df_fct[[vars$pt]] <- pt_2_factor(df_fct[[vars$pt]]) - return(df_fct) - - } + "as.data.frame", + signature = "DataLongitudinal", + definition = function(x, row.names, optional, ...) { + df_fct <- x@data + vars <- extractVariableNames(x) + df_fct[[vars$pt]] <- pt_2_factor(df_fct[[vars$pt]]) + return(df_fct) + } ) #' @rdname extractVariableNames setMethod( - f = "extractVariableNames", - signature = "DataLongitudinal", - definition = function(object) { - list( - pt = object@subject, - frm = object@formula, - time = as.character(object@formula[[3]]), - outcome = as.character(object@formula[[2]]), - threshold = object@threshold - ) - } + f = "extractVariableNames", + signature = "DataLongitudinal", + definition = function(object) { + list( + pt = object@subject, + frm = object@formula, + time = as.character(object@formula[[3]]), + outcome = as.character(object@formula[[2]]), + threshold = object@threshold + ) + } ) #' @export setMethod( - f = "as.list", - signature = "DataLongitudinal", - definition = function(x) { - - df <- as.data.frame(x) - vars <- extractVariableNames(x) - - mat_sld_index <- model.matrix( - as.formula(paste("~", vars$pt)), - data = df - ) |> - t() - - # TODO - Maybe reimplement this using a more robust approach than magic number - adj_threshold <- if (is.null(vars$threshold)) { - -999999 - } else { - vars$threshold - } - - index_obs <- which(df[[vars$outcome]] >= adj_threshold) - index_cen <- which(df[[vars$outcome]] < adj_threshold) - - sparse_mat_inds_all_y <- rstan::extract_sparse_parts(mat_sld_index) - sparse_mat_inds_obs_y <- rstan::extract_sparse_parts(mat_sld_index[, index_obs]) - sparse_mat_inds_cens_y <- rstan::extract_sparse_parts(mat_sld_index[, index_cen]) - - model_data <- list( - Nta_total = nrow(df), - Yobs = df[[vars$outcome]], - Tobs = df[[vars$time]], - Ythreshold = adj_threshold, - - # Number of individuals and tumor assessments. - Nta_obs_y = length(index_obs), - Nta_cens_y = length(index_cen), - - # Index vectors - ind_index = as.numeric(df[[vars$pt]]), - obs_y_index = index_obs, - cens_y_index = index_cen, - - # Sparse matrix parameters - # Matrix of individuals x observed tumor assessments. - n_mat_inds_obs_y = c( - length(sparse_mat_inds_obs_y$w), - length(sparse_mat_inds_obs_y$v), - length(sparse_mat_inds_obs_y$u) - ), - w_mat_inds_obs_y = sparse_mat_inds_obs_y$w, - v_mat_inds_obs_y = sparse_mat_inds_obs_y$v, - u_mat_inds_obs_y = sparse_mat_inds_obs_y$u, - - # Matrix of individuals x censored tumor assessments. - n_mat_inds_cens_y = c( - length(sparse_mat_inds_cens_y$w), - length(sparse_mat_inds_cens_y$v), - length(sparse_mat_inds_cens_y$u) - ), - w_mat_inds_cens_y = sparse_mat_inds_cens_y$w, - v_mat_inds_cens_y = sparse_mat_inds_cens_y$v, - u_mat_inds_cens_y = sparse_mat_inds_cens_y$u, - - # Matrix of all individuals tumour assessments - n_mat_inds_all_y = c( - length(sparse_mat_inds_all_y$w), - length(sparse_mat_inds_all_y$v), - length(sparse_mat_inds_all_y$u) - ), - w_mat_inds_all_y = sparse_mat_inds_all_y$w, - v_mat_inds_all_y = sparse_mat_inds_all_y$v, - u_mat_inds_all_y = sparse_mat_inds_all_y$u - ) - - return(model_data) + f = "as.list", + signature = "DataLongitudinal", + definition = function(x) { + df <- as.data.frame(x) + vars <- extractVariableNames(x) + + mat_sld_index <- model.matrix( + as.formula(paste("~", vars$pt)), + data = df + ) |> + t() + + # TODO - Maybe reimplement this using a more robust approach than magic number + adj_threshold <- if (is.null(vars$threshold)) { + -999999 + } else { + vars$threshold } -) + index_obs <- which(df[[vars$outcome]] >= adj_threshold) + index_cen <- which(df[[vars$outcome]] < adj_threshold) + + sparse_mat_inds_all_y <- rstan::extract_sparse_parts(mat_sld_index) + sparse_mat_inds_obs_y <- rstan::extract_sparse_parts(mat_sld_index[, index_obs]) + sparse_mat_inds_cens_y <- rstan::extract_sparse_parts(mat_sld_index[, index_cen]) + + model_data <- list( + Nta_total = nrow(df), + Yobs = df[[vars$outcome]], + Tobs = df[[vars$time]], + Ythreshold = adj_threshold, + + # Number of individuals and tumor assessments. + Nta_obs_y = length(index_obs), + Nta_cens_y = length(index_cen), + + # Index vectors + ind_index = as.numeric(df[[vars$pt]]), + obs_y_index = index_obs, + cens_y_index = index_cen, + + # Sparse matrix parameters + # Matrix of individuals x observed tumor assessments. + n_mat_inds_obs_y = c( + length(sparse_mat_inds_obs_y$w), + length(sparse_mat_inds_obs_y$v), + length(sparse_mat_inds_obs_y$u) + ), + w_mat_inds_obs_y = sparse_mat_inds_obs_y$w, + v_mat_inds_obs_y = sparse_mat_inds_obs_y$v, + u_mat_inds_obs_y = sparse_mat_inds_obs_y$u, + + # Matrix of individuals x censored tumor assessments. + n_mat_inds_cens_y = c( + length(sparse_mat_inds_cens_y$w), + length(sparse_mat_inds_cens_y$v), + length(sparse_mat_inds_cens_y$u) + ), + w_mat_inds_cens_y = sparse_mat_inds_cens_y$w, + v_mat_inds_cens_y = sparse_mat_inds_cens_y$v, + u_mat_inds_cens_y = sparse_mat_inds_cens_y$u, + + # Matrix of all individuals tumour assessments + n_mat_inds_all_y = c( + length(sparse_mat_inds_all_y$w), + length(sparse_mat_inds_all_y$v), + length(sparse_mat_inds_all_y$u) + ), + w_mat_inds_all_y = sparse_mat_inds_all_y$w, + v_mat_inds_all_y = sparse_mat_inds_all_y$v, + u_mat_inds_all_y = sparse_mat_inds_all_y$u + ) + return(model_data) + } +) diff --git a/R/DataSurvival.R b/R/DataSurvival.R index 002f4cd9..4ca510db 100644 --- a/R/DataSurvival.R +++ b/R/DataSurvival.R @@ -1,17 +1,16 @@ - #' @include generics.R NULL #' @rdname DataSurvival .DataSurvival <- setClass( - Class = "DataSurvival", - representation = list( - data = "data.frame", - formula = "formula", - subject = "character", - arm = "character", - study = "character" - ) + Class = "DataSurvival", + representation = list( + data = "data.frame", + formula = "formula", + subject = "character", + arm = "character", + study = "character" + ) ) @@ -35,14 +34,14 @@ NULL #' @seealso [`DataJoint`], [`DataLongitudinal`] #' #' @export -DataSurvival <- function(data, formula, subject, arm , study) { - .DataSurvival( - data = remove_missing_rows(data, formula, c(subject, arm, study)), - formula = formula, - subject = subject, - arm = arm, - study = study - ) +DataSurvival <- function(data, formula, subject, arm, study) { + .DataSurvival( + data = remove_missing_rows(data, formula, c(subject, arm, study)), + formula = formula, + subject = subject, + arm = arm, + study = study + ) } @@ -50,124 +49,122 @@ DataSurvival <- function(data, formula, subject, arm , study) { setValidity( - "DataSurvival", - method = function(object) { - - dat <- object@data - x <- extractVariableNames(object) - - dnames <- names(dat) - - if (length(x$frm) != 3) { - return("`formula` should be a 2 sided formula") - } - if (as.character(x$frm[[2]][[1]]) != "Surv") { - return("The LHS of `formula` should be a call to survival::Surv()") - } - if (length(x$pt) != 1) { - return("`subject` should be a length 1 character vector") - } - if (length(x$arm) != 1) { - return("`arm` should be a length 1 character vector") - } - if (length(x$study) != 1) { - return("`study` should be a length 1 character vector") - } - for (v in c(x$pt, x$arm, x$study, x$time, x$event)) { - if (! v %in% dnames) { - return(sprintf("Variable %s is not in `data`", x$pt)) - } - } - if (length(dat[[x$pt]]) != length(unique(dat[[x$pt]]))) { - return("Only 1 survival observation per subject is supported") - } - return(TRUE) + "DataSurvival", + method = function(object) { + dat <- object@data + x <- extractVariableNames(object) + + dnames <- names(dat) + + if (length(x$frm) != 3) { + return("`formula` should be a 2 sided formula") + } + if (as.character(x$frm[[2]][[1]]) != "Surv") { + return("The LHS of `formula` should be a call to survival::Surv()") + } + if (length(x$pt) != 1) { + return("`subject` should be a length 1 character vector") + } + if (length(x$arm) != 1) { + return("`arm` should be a length 1 character vector") + } + if (length(x$study) != 1) { + return("`study` should be a length 1 character vector") } + for (v in c(x$pt, x$arm, x$study, x$time, x$event)) { + if (!v %in% dnames) { + return(sprintf("Variable %s is not in `data`", x$pt)) + } + } + if (length(dat[[x$pt]]) != length(unique(dat[[x$pt]]))) { + return("Only 1 survival observation per subject is supported") + } + return(TRUE) + } ) #' @rdname extractVariableNames #' @export setMethod( - "extractVariableNames", - signature = "DataSurvival", - definition = function(object, ...) { - list( - arm = object@arm, - study = object@study, - pt = object@subject, - frm = object@formula, - time = as.character(object@formula[[2]][[2]]), - event = as.character(object@formula[[2]][[3]]) - ) - } + "extractVariableNames", + signature = "DataSurvival", + definition = function(object, ...) { + list( + arm = object@arm, + study = object@study, + pt = object@subject, + frm = object@formula, + time = as.character(object@formula[[2]][[2]]), + event = as.character(object@formula[[2]][[3]]) + ) + } ) #' @export setMethod( - "as.data.frame", - signature = "DataSurvival", - definition = function(x, row.names, optional, ...) { - as(x, "data.frame") - } + "as.data.frame", + signature = "DataSurvival", + definition = function(x, row.names, optional, ...) { + as(x, "data.frame") + } ) #' @export setAs( - from = "DataSurvival", - to = "data.frame", - def = function(from) { - df <- from@data - vars <- extractVariableNames(from) - df_fct <- df - df_fct[[vars$arm]] <- factor(df_fct[[vars$arm]]) - df_fct[[vars$study]] <- factor(df_fct[[vars$study]]) - df_fct[[vars$pt]] <- pt_2_factor(df_fct[[vars$pt]]) - return(df_fct) - } + from = "DataSurvival", + to = "data.frame", + def = function(from) { + df <- from@data + vars <- extractVariableNames(from) + df_fct <- df + df_fct[[vars$arm]] <- factor(df_fct[[vars$arm]]) + df_fct[[vars$study]] <- factor(df_fct[[vars$study]]) + df_fct[[vars$pt]] <- pt_2_factor(df_fct[[vars$pt]]) + return(df_fct) + } ) #' @export setMethod( - "as.list", - signature = "DataSurvival", - definition = function(x, ...) { - - df <- as(x, "data.frame") - vars <- extractVariableNames(x) - - design_mat <- stats::model.matrix(vars$frm, data = df) - remove_index <- grep("(Intercept)", colnames(design_mat), fixed = TRUE) - design_mat <- design_mat[, -remove_index, drop = FALSE] - - # Parameters for efficient integration of hazard function -> survival function - gh_parameters <- statmod::gauss.quad(n = 15, kind = "legendre") - - model_data <- list( - Nind = nrow(df), - Nind_dead = sum(df[[vars$event]]), - dead_ind_index = which(df[[vars$event]] == 1), - Times = df[[vars$time]], - p_os_cov_design = ncol(design_mat), - os_cov_design = design_mat, - n_nodes = length(gh_parameters$nodes), - nodes = gh_parameters$nodes, - weights = gh_parameters$weights, - study_index = as.numeric(df[[vars$study]]), - arm_index = as.numeric(df[[vars$arm]]) - ) - - arm_data <- get_arm_study_data( - as.numeric(df[[vars$pt]]), - as.numeric(df[[vars$arm]]), - as.numeric(df[[vars$study]]) - ) - - - return(append(model_data, arm_data)) - } + "as.list", + signature = "DataSurvival", + definition = function(x, ...) { + df <- as(x, "data.frame") + vars <- extractVariableNames(x) + + design_mat <- stats::model.matrix(vars$frm, data = df) + remove_index <- grep("(Intercept)", colnames(design_mat), fixed = TRUE) + design_mat <- design_mat[, -remove_index, drop = FALSE] + + # Parameters for efficient integration of hazard function -> survival function + gh_parameters <- statmod::gauss.quad(n = 15, kind = "legendre") + + model_data <- list( + Nind = nrow(df), + Nind_dead = sum(df[[vars$event]]), + dead_ind_index = which(df[[vars$event]] == 1), + Times = df[[vars$time]], + p_os_cov_design = ncol(design_mat), + os_cov_design = design_mat, + n_nodes = length(gh_parameters$nodes), + nodes = gh_parameters$nodes, + weights = gh_parameters$weights, + study_index = as.numeric(df[[vars$study]]), + arm_index = as.numeric(df[[vars$arm]]) + ) + + arm_data <- get_arm_study_data( + as.numeric(df[[vars$pt]]), + as.numeric(df[[vars$arm]]), + as.numeric(df[[vars$study]]) + ) + + + return(append(model_data, arm_data)) + } ) @@ -184,42 +181,41 @@ setMethod( #' #' @keywords internal get_arm_study_data <- function(pt, arm, study) { - - assert_that( - length(pt) == length(arm), - length(arm) == length(study), - msg = "`pt`, `arm` and `study` must all have the same length" - ) - - result <- list() - - u_arm_study <- unique(data.frame(study = study, arm = arm)) - - # TODO - Assumption that arms are unique to studies - assert_that( - length(u_arm_study[["arm"]]) == length(unique(u_arm_study[["arm"]])), - length(unique(arm)) == nrow(u_arm_study), - msg = "Arms must be unique across Studies e.g. arm A can't be in Study-1 and Study-2" - ) - - u_arm_study_ord <- u_arm_study[order(u_arm_study[["arm"]]), ] - - n_per_arm <- vector("numeric", length = length(u_arm_study_ord[["arm"]])) - index_per_arm <- vector("numeric", length = 0) - - uniq_arm <- u_arm_study_ord[["arm"]] - for (i in seq_along(uniq_arm)) { - arm_selected <- uniq_arm[[i]] - index_per_arm <- c(index_per_arm, which(arm == arm_selected)) - n_per_arm[[i]] <- sum(arm == arm_selected) - } - - result[["index_per_arm"]] <- index_per_arm - result[["n_index_per_arm"]] <- n_per_arm - result[["n_arms"]] <- nrow(u_arm_study) - result[["arm_to_study_index"]] <- as.numeric(u_arm_study[["study"]]) - result[["n_studies"]] <- length(unique(u_arm_study[["study"]])) - return(result) + assert_that( + length(pt) == length(arm), + length(arm) == length(study), + msg = "`pt`, `arm` and `study` must all have the same length" + ) + + result <- list() + + u_arm_study <- unique(data.frame(study = study, arm = arm)) + + # TODO - Assumption that arms are unique to studies + assert_that( + length(u_arm_study[["arm"]]) == length(unique(u_arm_study[["arm"]])), + length(unique(arm)) == nrow(u_arm_study), + msg = "Arms must be unique across Studies e.g. arm A can't be in Study-1 and Study-2" + ) + + u_arm_study_ord <- u_arm_study[order(u_arm_study[["arm"]]), ] + + n_per_arm <- vector("numeric", length = length(u_arm_study_ord[["arm"]])) + index_per_arm <- vector("numeric", length = 0) + + uniq_arm <- u_arm_study_ord[["arm"]] + for (i in seq_along(uniq_arm)) { + arm_selected <- uniq_arm[[i]] + index_per_arm <- c(index_per_arm, which(arm == arm_selected)) + n_per_arm[[i]] <- sum(arm == arm_selected) + } + + result[["index_per_arm"]] <- index_per_arm + result[["n_index_per_arm"]] <- n_per_arm + result[["n_arms"]] <- nrow(u_arm_study) + result[["arm_to_study_index"]] <- as.numeric(u_arm_study[["study"]]) + result[["n_studies"]] <- length(unique(u_arm_study[["study"]])) + return(result) } @@ -235,11 +231,8 @@ get_arm_study_data <- function(pt, arm, study) { #' @returns a `factor` vector with the levels in alphabetical order #' @keywords internal pt_2_factor <- function(pt) { - pt_char <- as.character(pt) - pt_uniq <- unique(pt_char) - pt_uniq_ord <- pt_uniq[order(pt_uniq)] - factor(pt_char, levels = pt_uniq_ord) + pt_char <- as.character(pt) + pt_uniq <- unique(pt_char) + pt_uniq_ord <- pt_uniq[order(pt_uniq)] + factor(pt_char, levels = pt_uniq_ord) } - - - diff --git a/R/JointModel.R b/R/JointModel.R index c715e971..14949ebb 100755 --- a/R/JointModel.R +++ b/R/JointModel.R @@ -1,4 +1,3 @@ - #' @include StanModule.R #' @include generics.R #' @include ParameterList.R @@ -6,43 +5,42 @@ NULL .JointModel <- setClass( - Class = "JointModel", - slots = list( - stan = "StanModule", - parameters = "ParameterList" - ) + Class = "JointModel", + slots = list( + stan = "StanModule", + parameters = "ParameterList" + ) ) #' @export JointModel <- function(longitudinal = NULL, survival = NULL, link = NULL) { - - longitudinal_linked <- addLink(longitudinal, link) - - parameters <- merge( - getParameters(longitudinal_linked), - getParameters(survival) - ) - - base_model <- paste0(read_stan("base/base.stan"), collapse = "\n") - - stan_full <- jinjar::render( - .x = base_model, - longditudinal = add_missing_stan_blocks(as.list(longitudinal_linked)), - survival = add_missing_stan_blocks(as.list(survival)), - priors = as.list(parameters), - link_none = class(link)[[1]] == "LinkNone" | is.null(link) - ) - - full_plus_funs <- merge( - StanModule("base/functions.stan"), - StanModule(stan_full) - ) - - .JointModel( - stan = full_plus_funs, - parameters = parameters - ) + longitudinal_linked <- addLink(longitudinal, link) + + parameters <- merge( + getParameters(longitudinal_linked), + getParameters(survival) + ) + + base_model <- paste0(read_stan("base/base.stan"), collapse = "\n") + + stan_full <- jinjar::render( + .x = base_model, + longditudinal = add_missing_stan_blocks(as.list(longitudinal_linked)), + survival = add_missing_stan_blocks(as.list(survival)), + priors = as.list(parameters), + link_none = class(link)[[1]] == "LinkNone" | is.null(link) + ) + + full_plus_funs <- merge( + StanModule("base/functions.stan"), + StanModule(stan_full) + ) + + .JointModel( + stan = full_plus_funs, + parameters = parameters + ) } @@ -51,97 +49,82 @@ JointModel <- function(longitudinal = NULL, survival = NULL, link = NULL) { #' @param x A `JointModel` object #' @export setMethod( - f = "as.character", - signature = "JointModel", - definition = function(x) { - as.character(x@stan) - } + f = "as.character", + signature = "JointModel", + definition = function(x) { + as.character(x@stan) + } ) #' @export setMethod( - f = "write_stan", - signature = "JointModel", - definition = function(x, file_path) { - fi <- file(file_path, open = "w") - writeLines(as.character(x), fi) - close(fi) - } + f = "write_stan", + signature = "JointModel", + definition = function(x, file_path) { + fi <- file(file_path, open = "w") + writeLines(as.character(x), fi) + close(fi) + } ) setMethod( - f = "compileStanModel", - signature = "JointModel", - definition = function(object, exe_file) { - x <- compileStanModel(object@stan, exe_file) - invisible(x) - } + f = "compileStanModel", + signature = "JointModel", + definition = function(object, exe_file) { + x <- compileStanModel(object@stan, exe_file) + invisible(x) + } ) setMethod( - f = "sampleStanModel", - signature = "JointModel", - definition = function(object, data, ..., exe_file = NULL) { - - args <- list(...) - - if (is(data, "DataJoint")) { - args[["data"]] <- as.list(data) - } else if (is(data, "list")) { - args[["data"]] <- data - } else { - stop("`data` must either be a list or a DataJoint object") - } - - if (!"init" %in% names(args)) { - values_initial <- initialValues(object) - values_sizes <- size(object@parameters) - values_sizes_complete <- replace_with_lookup(values_sizes, args[["data"]]) - values_initial_expanded <- expand_initial_values(values_initial, values_sizes_complete) - args[["init"]] <- function() values_initial_expanded - } - - model <- compileStanModel(object, exe_file) - do.call(model$sample, args) + f = "sampleStanModel", + signature = "JointModel", + definition = function(object, data, ..., exe_file = NULL) { + args <- list(...) + + if (is(data, "DataJoint")) { + args[["data"]] <- as.list(data) + } else if (is(data, "list")) { + args[["data"]] <- data + } else { + stop("`data` must either be a list or a DataJoint object") } + + if (!"init" %in% names(args)) { + values_initial <- initialValues(object) + values_sizes <- size(object@parameters) + values_sizes_complete <- replace_with_lookup(values_sizes, args[["data"]]) + values_initial_expanded <- expand_initial_values(values_initial, values_sizes_complete) + args[["init"]] <- function() values_initial_expanded + } + + model <- compileStanModel(object, exe_file) + do.call(model$sample, args) + } ) #' @export setMethod( - f = "initialValues", - signature = "JointModel", - definition = function(object) { - initialValues(object@parameters) - } + f = "initialValues", + signature = "JointModel", + definition = function(object) { + initialValues(object@parameters) + } ) add_missing_stan_blocks <- function(x) { - # STAN_BLOCKS is defined as a global variable in stan_module.R - # TODO - Make it an argument to the function - for (block in names(STAN_BLOCKS)) { - if (is.null(x[[block]])) { - x[[block]] <- "" - } + # STAN_BLOCKS is defined as a global variable in stan_module.R + # TODO - Make it an argument to the function + for (block in names(STAN_BLOCKS)) { + if (is.null(x[[block]])) { + x[[block]] <- "" } - return(x) + } + return(x) } - - - - - - - - - - - - - - diff --git a/R/Link.R b/R/Link.R index 85887b00..bc00a9eb 100755 --- a/R/Link.R +++ b/R/Link.R @@ -4,52 +4,48 @@ NULL .Link <- setClass( - Class = "Link", - slots = list( - "stan" = "StanModule", - "parameters" = "ParameterList" - ) + Class = "Link", + slots = list( + "stan" = "StanModule", + "parameters" = "ParameterList" + ) ) #' @export Link <- function(stan = StanModule(), parameters = ParameterList(), ...) { - .Link(stan = stan, parameters = parameters, ...) + .Link(stan = stan, parameters = parameters, ...) } setMethod( - f = "addLink", - signature = c("LongitudinalModel", "Link"), - definition = function(x, y, ...) { - x@stan <- merge(x@stan, y@stan) - x@parameters <- merge(x@parameters, y@parameters) - x - } + f = "addLink", + signature = c("LongitudinalModel", "Link"), + definition = function(x, y, ...) { + x@stan <- merge(x@stan, y@stan) + x@parameters <- merge(x@parameters, y@parameters) + x + } ) #' @export setMethod( - f = "initialValues", - signature = "Link", - definition = function(object) { - initialValues(object@parameters) - } + f = "initialValues", + signature = "Link", + definition = function(object) { + initialValues(object@parameters) + } ) -#' @export +#' @export setMethod( - f = "as.StanModule", - signature = "Link", - definition = function(object) { - object@stan - } + f = "as.StanModule", + signature = "Link", + definition = function(object) { + object@stan + } ) - - - - diff --git a/R/LinkNone.R b/R/LinkNone.R index 00302dca..55802ddf 100755 --- a/R/LinkNone.R +++ b/R/LinkNone.R @@ -1,4 +1,3 @@ - #' @include generics.R #' @include Link.R #' @include LongitudinalModel.R @@ -6,22 +5,21 @@ NULL .LinkNone <- setClass( - Class = "LinkNone", - contains = "Link" + Class = "LinkNone", + contains = "Link" ) #' @export LinkNone <- function() { - stan = StanModule() - .LinkNone(Link(stan = stan)) + stan <- StanModule() + .LinkNone(Link(stan = stan)) } setMethod( - f = "addLink", - signature = c("LongitudinalModel", "LinkNone"), - definition = function(x, y, ...) { - x - } + f = "addLink", + signature = c("LongitudinalModel", "LinkNone"), + definition = function(x, y, ...) { + x + } ) - diff --git a/R/LongitudinalGSF.R b/R/LongitudinalGSF.R index 1a6ccf58..95e69043 100755 --- a/R/LongitudinalGSF.R +++ b/R/LongitudinalGSF.R @@ -1,4 +1,3 @@ - #' @include LongitudinalModel.R #' @include StanModule.R #' @include generics.R @@ -8,13 +7,13 @@ NULL .LongitudinalGSF <- setClass( - Class = "LongitudinalGSF", - contains = "LongitudinalModel" + Class = "LongitudinalGSF", + contains = "LongitudinalModel" ) #' @export -LongitudinalGSF <- function ( +LongitudinalGSF <- function( mu_bsld = prior_lognormal(log(55), 5, init = 55), mu_ks = prior_lognormal(0, 0.5, init = 0), mu_kg = prior_lognormal(-0.36, 1, init = -0.36), @@ -27,140 +26,133 @@ LongitudinalGSF <- function ( tilde_bsld = prior_normal(0, 5, init = 0), tilde_ks = prior_normal(0, 5, init = 0), tilde_kg = prior_normal(0, 5, init = 0), - tilde_phi = prior_normal(0, 5, init = 0) -) { - x <- LongitudinalModel( - stan = merge( - StanModule("lm-gsf/model.stan"), - StanModule("lm-gsf/functions.stan") - ), - parameters = ParameterList( - Parameter(name = "lm_gsf_mu_bsld", prior = mu_bsld, size = "n_studies"), - Parameter(name = "lm_gsf_mu_ks", prior = mu_ks, size = "n_arms"), - Parameter(name = "lm_gsf_mu_kg", prior = mu_kg, size = "n_arms"), - Parameter(name = "lm_gsf_mu_phi", prior = mu_phi, size = "n_arms"), - Parameter(name = "lm_gsf_omega_bsld", prior = omega_bsld, size = 1), - Parameter(name = "lm_gsf_omega_ks", prior = omega_ks, size = 1), - Parameter(name = "lm_gsf_omega_kg", prior = omega_kg, size = 1), - Parameter(name = "lm_gsf_omega_phi", prior = omega_phi, size = 1), - Parameter(name = "lm_gsf_sigma", prior = sigma, size = 1), - Parameter(name = "lm_gsf_eta_tilde_bsld", prior = tilde_bsld, size = "Nind"), - Parameter(name = "lm_gsf_eta_tilde_ks", prior = tilde_ks, size = "Nind"), - Parameter(name = "lm_gsf_eta_tilde_kg", prior = tilde_kg, size = "Nind"), - Parameter(name = "lm_gsf_eta_tilde_phi", prior = tilde_phi, size = "Nind") - ) + tilde_phi = prior_normal(0, 5, init = 0)) { + x <- LongitudinalModel( + stan = merge( + StanModule("lm-gsf/model.stan"), + StanModule("lm-gsf/functions.stan") + ), + parameters = ParameterList( + Parameter(name = "lm_gsf_mu_bsld", prior = mu_bsld, size = "n_studies"), + Parameter(name = "lm_gsf_mu_ks", prior = mu_ks, size = "n_arms"), + Parameter(name = "lm_gsf_mu_kg", prior = mu_kg, size = "n_arms"), + Parameter(name = "lm_gsf_mu_phi", prior = mu_phi, size = "n_arms"), + Parameter(name = "lm_gsf_omega_bsld", prior = omega_bsld, size = 1), + Parameter(name = "lm_gsf_omega_ks", prior = omega_ks, size = 1), + Parameter(name = "lm_gsf_omega_kg", prior = omega_kg, size = 1), + Parameter(name = "lm_gsf_omega_phi", prior = omega_phi, size = 1), + Parameter(name = "lm_gsf_sigma", prior = sigma, size = 1), + Parameter(name = "lm_gsf_eta_tilde_bsld", prior = tilde_bsld, size = "Nind"), + Parameter(name = "lm_gsf_eta_tilde_ks", prior = tilde_ks, size = "Nind"), + Parameter(name = "lm_gsf_eta_tilde_kg", prior = tilde_kg, size = "Nind"), + Parameter(name = "lm_gsf_eta_tilde_phi", prior = tilde_phi, size = "Nind") ) - .LongitudinalGSF(x) + ) + .LongitudinalGSF(x) } .LinkGSF <- setClass( - Class = "LinkGSF", - contains = "Link" + Class = "LinkGSF", + contains = "Link" ) #' @export LinkGSF <- function( components = list( - link_gsf_dsld(), - link_gsf_ttg() - ) -) { - - items <- lapply( - components, - function(x) { - list( - parameter = x@parameter_name, - contribution_function = x@contribution_fname - ) - } - ) - - rendered_link <- jinjar::render( - .x = paste0(read_stan("lm-gsf/link.stan"), collapse = "\n"), - items = items - ) - - parameters <- ParameterList() - stan_components <- StanModule() - for (item in components) { - parameters <- merge(parameters, item@parameter) - stan_components <- merge(stan_components, item@stan) + link_gsf_dsld(), + link_gsf_ttg() + )) { + items <- lapply( + components, + function(x) { + list( + parameter = x@parameter_name, + contribution_function = x@contribution_fname + ) } - - stan_full <- merge( - stan_components, - StanModule(rendered_link) - ) - - x <- Link( - stan = stan_full, - parameters = parameters - ) - - .LinkGSF(x) + ) + + rendered_link <- jinjar::render( + .x = paste0(read_stan("lm-gsf/link.stan"), collapse = "\n"), + items = items + ) + + parameters <- ParameterList() + stan_components <- StanModule() + for (item in components) { + parameters <- merge(parameters, item@parameter) + stan_components <- merge(stan_components, item@stan) + } + + stan_full <- merge( + stan_components, + StanModule(rendered_link) + ) + + x <- Link( + stan = stan_full, + parameters = parameters + ) + + .LinkGSF(x) } .link_gsf_abstract <- setClass( - Class = "link_gsf_abstract", - slots = list( - "stan" = "StanModule", - "parameter" = "ParameterList", - "parameter_name" = "character", - "contribution_fname" = "character" - ) + Class = "link_gsf_abstract", + slots = list( + "stan" = "StanModule", + "parameter" = "ParameterList", + "parameter_name" = "character", + "contribution_fname" = "character" + ) ) #' @export link_gsf_abstract <- function(stan, parameter, parameter_name, contribution_fname) { - .link_gsf_abstract( - parameter = parameter, - parameter_name = parameter_name, - contribution_fname = contribution_fname, - stan = StanModule( - paste0( - "functions {\n", - as.list(stan)[["functions"]], - "\n}" - ) - ) + .link_gsf_abstract( + parameter = parameter, + parameter_name = parameter_name, + contribution_fname = contribution_fname, + stan = StanModule( + paste0( + "functions {\n", + as.list(stan)[["functions"]], + "\n}" + ) ) + ) } .link_gsf_ttg <- setClass( - Class = "link_gsf_ttg", - contains = "link_gsf_abstract" + Class = "link_gsf_ttg", + contains = "link_gsf_abstract" ) #' @export link_gsf_ttg <- function( - gamma = prior_normal(0, 5, init = 0) -) { - link_gsf_abstract( - stan = StanModule("lm-gsf/link_ttg.stan"), - parameter = ParameterList(Parameter(name = "lm_gsf_gamma", prior = gamma, size = 1)), - parameter_name = "lm_gsf_gamma", - contribution_fname = "link_ttg_contribution" - ) + gamma = prior_normal(0, 5, init = 0)) { + link_gsf_abstract( + stan = StanModule("lm-gsf/link_ttg.stan"), + parameter = ParameterList(Parameter(name = "lm_gsf_gamma", prior = gamma, size = 1)), + parameter_name = "lm_gsf_gamma", + contribution_fname = "link_ttg_contribution" + ) } .link_gsf_dsld <- setClass( - Class = "link_gsf_dsld", - contains = "link_gsf_abstract" + Class = "link_gsf_dsld", + contains = "link_gsf_abstract" ) #' @export link_gsf_dsld <- function( - beta = Parameter(prior_normal(0, 5), init = 0) -) { - link_gsf_abstract( - stan = StanModule("lm-gsf/link_dsld.stan"), - parameter = ParameterList(Parameter(name = "lm_gsf_beta", prior = beta, size = 1)), - parameter_name = "lm_gsf_beta", - contribution_fname = "link_dsld_contribution" - ) + beta = Parameter(prior_normal(0, 5), init = 0)) { + link_gsf_abstract( + stan = StanModule("lm-gsf/link_dsld.stan"), + parameter = ParameterList(Parameter(name = "lm_gsf_beta", prior = beta, size = 1)), + parameter_name = "lm_gsf_beta", + contribution_fname = "link_dsld_contribution" + ) } - - diff --git a/R/LongitudinalModel.R b/R/LongitudinalModel.R index a30b978e..4462484a 100755 --- a/R/LongitudinalModel.R +++ b/R/LongitudinalModel.R @@ -1,30 +1,23 @@ - #' @include StanModel.R NULL .LongitudinalModel <- setClass( - Class = "LongitudinalModel", - contains = "StanModel" + Class = "LongitudinalModel", + contains = "StanModel" ) #' @export LongitudinalModel <- function(stan = StanModule(), parameters = ParameterList(), ...) { - - base_long <- StanModule( - x = "base/longitudinal.stan" - ) - - .LongitudinalModel( - StanModel( - stan = merge(base_long, stan), - parameters = parameters, - ... - ) + base_long <- StanModule( + x = "base/longitudinal.stan" + ) + + .LongitudinalModel( + StanModel( + stan = merge(base_long, stan), + parameters = parameters, + ... ) + ) } - - - - - diff --git a/R/LongitudinalRandomSlope.R b/R/LongitudinalRandomSlope.R index 47d50156..4b178bad 100755 --- a/R/LongitudinalRandomSlope.R +++ b/R/LongitudinalRandomSlope.R @@ -1,12 +1,11 @@ - #' @include LongitudinalModel.R #' @include Link.R NULL .LongitudinalRandomSlope <- setClass( - Class = "LongitudinalRandomSlope", - contains = "LongitudinalModel" + Class = "LongitudinalRandomSlope", + contains = "LongitudinalModel" ) @@ -16,55 +15,47 @@ LongitudinalRandomSlope <- function( slope_mu = prior_normal(0, 15, init = 0.001), slope_sigma = prior_cauchy(0, 4, init = 0.4), sigma = prior_cauchy(0, 4, init = 0.4), - random_slope = prior_none(init = 0) -) { - - stan <- StanModule( - x = "lm-random-slope/model.stan" - ) - - if (!as.character(random_slope) == "") { - stop("`random_slope` must be a `prior_none()`") - } - - .LongitudinalRandomSlope( - LongitudinalModel( - stan = stan, - parameters = ParameterList( - Parameter(name = "lm_rs_intercept", prior = intercept, size = 1), - Parameter(name = "lm_rs_slope_mu", prior = slope_mu, size = "n_arms"), - Parameter(name = "lm_rs_slope_sigma", prior = slope_sigma, size = 1), - Parameter(name = "lm_rs_sigma", prior = sigma, size = 1), - Parameter(name = "lm_rs_ind_rnd_slope", prior = random_slope, size = "Nind") - ) - ) + random_slope = prior_none(init = 0)) { + stan <- StanModule( + x = "lm-random-slope/model.stan" + ) + + if (!as.character(random_slope) == "") { + stop("`random_slope` must be a `prior_none()`") + } + + .LongitudinalRandomSlope( + LongitudinalModel( + stan = stan, + parameters = ParameterList( + Parameter(name = "lm_rs_intercept", prior = intercept, size = 1), + Parameter(name = "lm_rs_slope_mu", prior = slope_mu, size = "n_arms"), + Parameter(name = "lm_rs_slope_sigma", prior = slope_sigma, size = 1), + Parameter(name = "lm_rs_sigma", prior = sigma, size = 1), + Parameter(name = "lm_rs_ind_rnd_slope", prior = random_slope, size = "Nind") + ) ) + ) } .LinkRandomSlope <- setClass( - Class = "LinkRandomSlope", - contains = "Link" + Class = "LinkRandomSlope", + contains = "Link" ) #' @export LinkRandomSlope <- function( - link_lm_phi = prior_normal(0.2, 0.5, init = 0.02) -) { - .LinkRandomSlope( - Link( - stan = StanModule( - x = "lm-random-slope/links.stan" - ), - parameters = ParameterList( - Parameter(name = "link_lm_phi", prior = link_lm_phi, size = 1) - ) - ) + link_lm_phi = prior_normal(0.2, 0.5, init = 0.02)) { + .LinkRandomSlope( + Link( + stan = StanModule( + x = "lm-random-slope/links.stan" + ), + parameters = ParameterList( + Parameter(name = "link_lm_phi", prior = link_lm_phi, size = 1) + ) ) + ) } - - - - - diff --git a/R/Parameter.R b/R/Parameter.R index 07044dfe..bc64a4fe 100755 --- a/R/Parameter.R +++ b/R/Parameter.R @@ -1,4 +1,3 @@ - #' @include generics.R #' @include Prior.R NULL @@ -7,87 +6,86 @@ NULL setClassUnion(name = "numeric_OR_character", c("numeric", "character")) .Parameter <- setClass( - Class = "Parameter", - slots = list( - "name" = "character", - "prior" = "Prior", - "size" = "numeric_OR_character" - ) + Class = "Parameter", + slots = list( + "name" = "character", + "prior" = "Prior", + "size" = "numeric_OR_character" + ) ) #' @export Parameter <- function(prior, name, size = 1) { - .Parameter( - prior = prior, - name = name, - size = size - ) + .Parameter( + prior = prior, + name = name, + size = size + ) } setValidity( - Class = "Parameter", - method = function(object) { - if (!length(object@name) == 1) { - return("Name must be a length 1 character vector") - } - if (is.character(object@size)) { - if (!length(object@size) == 1) { - return("Size must be a numeric vector or length 1 character vector") - } - } - return(TRUE) + Class = "Parameter", + method = function(object) { + if (!length(object@name) == 1) { + return("Name must be a length 1 character vector") + } + if (is.character(object@size)) { + if (!length(object@size) == 1) { + return("Size must be a numeric vector or length 1 character vector") + } } + return(TRUE) + } ) #' @export setMethod( - f = "as.character", - signature = "Parameter", - definition = function(x) as(x, "character") + f = "as.character", + signature = "Parameter", + definition = function(x) as(x, "character") ) #' @export setAs( - from = "Parameter", - to = "character", - def = function(from) { - if (as.character(from@prior) =="") { - return("") - } - glue::glue("{name} ~ {dist}", name = from@name, dist = as(from@prior, "character")) + from = "Parameter", + to = "character", + def = function(from) { + if (as.character(from@prior) == "") { + return("") } + glue::glue("{name} ~ {dist}", name = from@name, dist = as(from@prior, "character")) + } ) -#' @export +#' @export setMethod( - f = "names", - signature = "Parameter", - definition = function(x) x@name + f = "names", + signature = "Parameter", + definition = function(x) x@name ) -#' @export +#' @export setMethod( - f = "initialValues", - signature = "Parameter", - definition = function(object) initialValues(object@prior) + f = "initialValues", + signature = "Parameter", + definition = function(object) initialValues(object@prior) ) #' @export setMethod( - f = "size", - signature = "Parameter", - definition = function(object) object@size + f = "size", + signature = "Parameter", + definition = function(object) object@size ) - diff --git a/R/ParameterList.R b/R/ParameterList.R index b7a9c3aa..ac58c922 100644 --- a/R/ParameterList.R +++ b/R/ParameterList.R @@ -7,132 +7,131 @@ NULL .ParameterList <- setClass( - Class = "ParameterList", - slots = c( - parameters = "list" - ) + Class = "ParameterList", + slots = c( + parameters = "list" + ) ) #' @export ParameterList <- function(...) { - .ParameterList(parameters = list(...)) + .ParameterList(parameters = list(...)) } #' @export setValidity( - Class = "ParameterList", - method = function(object) { - is_parameters <- vapply(object@parameters, function(x) is(x, "Parameter"), logical(1)) - if (!all(is_parameters)) { - return("all elements must be of class 'Parameter'") - } - return(TRUE) + Class = "ParameterList", + method = function(object) { + is_parameters <- vapply(object@parameters, function(x) is(x, "Parameter"), logical(1)) + if (!all(is_parameters)) { + return("all elements must be of class 'Parameter'") } + return(TRUE) + } ) #' @export setAs( - from = "ParameterList", - to = "character", - def = function(from) { - strings <- vapply( - from@parameters, - as.character, - character(1) - ) - indentation <- paste0(rep(" ", 4), collapse = "") - strings_indented <- paste0(indentation, strings) - paste(strings_indented, collapse = "\n") - } + from = "ParameterList", + to = "character", + def = function(from) { + strings <- vapply( + from@parameters, + as.character, + character(1) + ) + indentation <- paste0(rep(" ", 4), collapse = "") + strings_indented <- paste0(indentation, strings) + paste(strings_indented, collapse = "\n") + } ) #' @export setMethod( - f = "as.character", - signature = "ParameterList", - definition = function(x) { - as(x, "character") - } + f = "as.character", + signature = "ParameterList", + definition = function(x) { + as(x, "character") + } ) #' @export setMethod( - f = "as.StanModule", - signature = "ParameterList", - definition = function(object) { - x <- paste( - "model {", - as(object, "character"), - "}", - sep = "\n" - ) - StanModule(x = x) - } + f = "as.StanModule", + signature = "ParameterList", + definition = function(object) { + x <- paste( + "model {", + as(object, "character"), + "}", + sep = "\n" + ) + StanModule(x = x) + } ) #' @export setMethod( - f = "merge", - signature = c(x = "ParameterList", y = "ParameterList"), - definition = function(x, y) { - parameters <- append(x@parameters, y@parameters) - do.call(ParameterList, parameters) - } + f = "merge", + signature = c(x = "ParameterList", y = "ParameterList"), + definition = function(x, y) { + parameters <- append(x@parameters, y@parameters) + do.call(ParameterList, parameters) + } ) #' @export setMethod( - f = "as.list", - signature = "ParameterList", - definition = function(x) { - as.list(as.StanModule(x)) - } + f = "as.list", + signature = "ParameterList", + definition = function(x) { + as.list(as.StanModule(x)) + } ) #' @export setMethod( - f = "initialValues", - signature = "ParameterList", - definition = function(object) { - vals <- lapply(object@parameters, initialValues) - name <- vapply(object@parameters, names, character(1)) - names(vals) <- name - return(vals) - } + f = "initialValues", + signature = "ParameterList", + definition = function(object) { + vals <- lapply(object@parameters, initialValues) + name <- vapply(object@parameters, names, character(1)) + names(vals) <- name + return(vals) + } ) -#' @export +#' @export setMethod( - f = "names", - signature = "ParameterList", - definition = function(x) { - vapply(x@parameters, names, character(1)) - } + f = "names", + signature = "ParameterList", + definition = function(x) { + vapply(x@parameters, names, character(1)) + } ) #' @export setMethod( - f = "size", - signature = "ParameterList", - definition = function(object) { - x <- lapply(object@parameters, size) - names(x) <- names(object) - return(x) - } + f = "size", + signature = "ParameterList", + definition = function(object) { + x <- lapply(object@parameters, size) + names(x) <- names(object) + return(x) + } ) - diff --git a/R/Prior.R b/R/Prior.R index 21aef8f7..11a171d8 100755 --- a/R/Prior.R +++ b/R/Prior.R @@ -1,99 +1,97 @@ - #' @include generics.R -NULL +NULL .Prior <- setClass( - Class = "Prior", - slots = c( - "parameters" = "list", - "repr" = "character", - "init" = "numeric" - ) + Class = "Prior", + slots = c( + "parameters" = "list", + "repr" = "character", + "init" = "numeric" + ) ) #' @export setMethod( - f = "as.character", - signature = "Prior", - definition = function(x) { - as(x, "character") - } + f = "as.character", + signature = "Prior", + definition = function(x) { + as(x, "character") + } ) #' @export setAs( - from = "Prior", - to = "character", - def = function(from) { - glue::glue(from@repr, .envir = list2env(from@parameters)) - } + from = "Prior", + to = "character", + def = function(from) { + glue::glue(from@repr, .envir = list2env(from@parameters)) + } ) #' @export prior_normal <- function(mu, sigma, init = mu) { - .Prior( - parameters = list(mu = mu, sigma = sigma), - repr = "normal({mu}, {sigma});", - init = init - ) + .Prior( + parameters = list(mu = mu, sigma = sigma), + repr = "normal({mu}, {sigma});", + init = init + ) } #' @export prior_cauchy <- function(mu, sigma, init = mu) { - .Prior( - parameters = list(mu = mu, sigma = sigma), - repr = "cauchy({mu}, {sigma});", - init = init - ) + .Prior( + parameters = list(mu = mu, sigma = sigma), + repr = "cauchy({mu}, {sigma});", + init = init + ) } #' @export -prior_gamma <- function(alpha, beta, init = alpha/beta) { - .Prior( - parameters = list(alpha = alpha, beta = beta), - repr = "gamma({alpha}, {beta});", - init = init - ) +prior_gamma <- function(alpha, beta, init = alpha / beta) { + .Prior( + parameters = list(alpha = alpha, beta = beta), + repr = "gamma({alpha}, {beta});", + init = init + ) } #' @export -prior_lognormal <- function(mu, sigma, init = exp(mu + (sigma^2)/2)) { - .Prior( - parameters = list(mu = mu, sigma = sigma), - repr = "lognormal({mu}, {sigma});", - init = init - ) +prior_lognormal <- function(mu, sigma, init = exp(mu + (sigma^2) / 2)) { + .Prior( + parameters = list(mu = mu, sigma = sigma), + repr = "lognormal({mu}, {sigma});", + init = init + ) } #' @export -prior_beta <- function(a, b, init = a/(a+b)) { - .Prior( - parameters = list(a = a, b = b), - repr = "beta({a}, {b});", - init = init - ) +prior_beta <- function(a, b, init = a / (a + b)) { + .Prior( + parameters = list(a = a, b = b), + repr = "beta({a}, {b});", + init = init + ) } #' @export prior_none <- function(init = 0.00001) { - .Prior( - parameters = list(), - repr = "", - init = init - ) + .Prior( + parameters = list(), + repr = "", + init = init + ) } -#' @export +#' @export setMethod( - f = "initialValues", - signature = "Prior", - definition = function(object) object@init + f = "initialValues", + signature = "Prior", + definition = function(object) object@init ) - diff --git a/R/StanModel.R b/R/StanModel.R index bb2c7bd1..c40e2751 100644 --- a/R/StanModel.R +++ b/R/StanModel.R @@ -1,42 +1,39 @@ - #' @include StanModule.R #' @include ParameterList.R NULL .StanModel <- setClass( - Class = "StanModel", - slots = list( - "stan" = "StanModule", - "parameters" = "ParameterList" - ) + Class = "StanModel", + slots = list( + "stan" = "StanModule", + "parameters" = "ParameterList" + ) ) StanModel <- function(stan, parameters, ...) { - .StanModel( - stan = stan, - parameters = parameters, - ... - ) + .StanModel( + stan = stan, + parameters = parameters, + ... + ) } setMethod( - f = "getParameters", - signature = "StanModel", - definition = function(object) { - object@parameters - } + f = "getParameters", + signature = "StanModel", + definition = function(object) { + object@parameters + } ) #' @export setMethod( - f = "as.list", - signature = c("StanModel"), - definition = function(x) { - as.list(x@stan) - } + f = "as.list", + signature = c("StanModel"), + definition = function(x) { + as.list(x@stan) + } ) - - diff --git a/R/StanModule.R b/R/StanModule.R index e1eefe95..d6a5b3b8 100755 --- a/R/StanModule.R +++ b/R/StanModule.R @@ -1,15 +1,14 @@ - #' @include generics.R NULL STAN_BLOCKS <- list( - functions = "functions", - data = "data", - transformed_data = "transformed data", - parameters = "parameters", - transformed_parameters = "transformed parameters", - model = "model", - generated_quantities = "generated quantities" + functions = "functions", + data = "data", + transformed_data = "transformed data", + parameters = "parameters", + transformed_parameters = "transformed parameters", + model = "model", + generated_quantities = "generated quantities" ) @@ -18,18 +17,18 @@ STAN_BLOCKS <- list( #' StanAll class object .StanModule <- setClass( - Class = "StanModule", - slots = list( - functions = "character", - data = "character", - transformed_data = "character", - parameters = "character", - transformed_parameters = "character", - model = "character", - generated_quantities = "character", - priors = "list", - inits = "list" - ) + Class = "StanModule", + slots = list( + functions = "character", + data = "character", + transformed_data = "character", + parameters = "character", + transformed_parameters = "character", + model = "character", + generated_quantities = "character", + priors = "list", + inits = "list" + ) ) @@ -38,27 +37,26 @@ StanModule <- function( x = "", priors = list(), inits = list(), + ...) { + assert_that( + is.character(x), + length(x) == 1, + msg = "`x` must be either a length 1 character vector" + ) + code <- read_stan(x) + code_fragments <- as_stan_fragments(code) + .StanModule( + functions = code_fragments$functions, + data = code_fragments$data, + transformed_data = code_fragments$transformed_data, + parameters = code_fragments$parameters, + transformed_parameters = code_fragments$transformed_parameters, + model = code_fragments$model, + generated_quantities = code_fragments$generated_quantities, + priors = priors, + inits = inits, ... -) { - assert_that( - is.character(x), - length(x) == 1, - msg = "`x` must be either a length 1 character vector" - ) - code <- read_stan(x) - code_fragments <- as_stan_fragments(code) - .StanModule( - functions = code_fragments$functions, - data = code_fragments$data, - transformed_data = code_fragments$transformed_data, - parameters = code_fragments$parameters, - transformed_parameters = code_fragments$transformed_parameters, - model = code_fragments$model, - generated_quantities = code_fragments$generated_quantities, - priors = priors, - inits = inits, - ... - ) + ) } @@ -92,73 +90,73 @@ StanModule <- function( #' @param x A `StanModule` object #' @export setMethod( - f = "as.character", - signature = "StanModule", - definition = function(x) { - as_stan_file( - functions = x@functions, - data = x@data, - transformed_data = x@transformed_data, - parameters = x@parameters, - transformed_parameters = x@transformed_parameters, - model = x@model, - generated_quantities = x@generated_quantities - ) - } + f = "as.character", + signature = "StanModule", + definition = function(x) { + as_stan_file( + functions = x@functions, + data = x@data, + transformed_data = x@transformed_data, + parameters = x@parameters, + transformed_parameters = x@transformed_parameters, + model = x@model, + generated_quantities = x@generated_quantities + ) + } ) #' @export setMethod( - f = "merge", - signature = c("StanModule", "StanModule"), - definition = function(x, y, ...) { - stan_fragments <- lapply( - names(STAN_BLOCKS), - \(par) { - new_string <- c(slot(x, par), slot(y, par)) - if (all(new_string == "")) { - return("") - } - return(new_string) - } - ) - names(stan_fragments) <- names(STAN_BLOCKS) - stan_code <- do.call(as_stan_file, stan_fragments) - StanModule( - x = stan_code, - priors = append(x@priors, y@priors), - inits = append(x@inits, y@inits) - ) - } + f = "merge", + signature = c("StanModule", "StanModule"), + definition = function(x, y, ...) { + stan_fragments <- lapply( + names(STAN_BLOCKS), + \(par) { + new_string <- c(slot(x, par), slot(y, par)) + if (all(new_string == "")) { + return("") + } + return(new_string) + } + ) + names(stan_fragments) <- names(STAN_BLOCKS) + stan_code <- do.call(as_stan_file, stan_fragments) + StanModule( + x = stan_code, + priors = append(x@priors, y@priors), + inits = append(x@inits, y@inits) + ) + } ) #' @export setMethod( - f = "compileStanModel", - signature = signature(object = "StanModule", exe_file = "character"), - definition = function(object, exe_file) { - if (!dir.exists(dirname(exe_file))) { - dir.create(dirname(exe_file), recursive = TRUE) - } - x <- cmdstanr::cmdstan_model( - stan_file = cmdstanr::write_stan_file(as.character(object)), - exe_file = exe_file - ) - invisible(x) + f = "compileStanModel", + signature = signature(object = "StanModule", exe_file = "character"), + definition = function(object, exe_file) { + if (!dir.exists(dirname(exe_file))) { + dir.create(dirname(exe_file), recursive = TRUE) } + x <- cmdstanr::cmdstan_model( + stan_file = cmdstanr::write_stan_file(as.character(object)), + exe_file = exe_file + ) + invisible(x) + } ) #' @export setMethod( - f = "compileStanModel", - signature = signature(object = "StanModule", exe_file = "empty"), - definition = function(object, exe_file) { - exe_file <- file.path(tempdir(), "model") - invisible(compileStanModel(object, exe_file)) - } + f = "compileStanModel", + signature = signature(object = "StanModule", exe_file = "empty"), + definition = function(object, exe_file) { + exe_file <- file.path(tempdir(), "model") + invisible(compileStanModel(object, exe_file)) + } ) @@ -167,11 +165,11 @@ setMethod( #' @export setMethod( - f = "show", - signature = "StanModule", - definition = function(object) { - print("StanModule Object") - } + f = "show", + signature = "StanModule", + definition = function(object) { + print("StanModule Object") + } ) @@ -179,16 +177,16 @@ setMethod( #' @export setMethod( - f = "as.list", - signature = "StanModule", - definition = function(x) { - string <- as.character(x) - li <- as_stan_fragments(string) - for (block in names(STAN_BLOCKS)) { - li[[block]] <- paste(li[[block]], collapse = "\n") - } - return(li) + f = "as.list", + signature = "StanModule", + definition = function(x) { + string <- as.character(x) + li <- as_stan_fragments(string) + for (block in names(STAN_BLOCKS)) { + li[[block]] <- paste(li[[block]], collapse = "\n") } + return(li) + } ) @@ -201,21 +199,21 @@ setMethod( #' @param filename A character string #' @importFrom assertthat assert_that is_file <- function(filename = NULL) { - if (is.null(filename)) { - return(FALSE) - } - assert_that( - is.character(filename), - length(filename) == 1, - msg = "`filename` must be a length 1 character" - ) - if (nchar(filename) > 1000) { - return(FALSE) - } - if (is.na(filename)) { - return(FALSE) - } - return(file.exists(filename) & !dir.exists(filename)) + if (is.null(filename)) { + return(FALSE) + } + assert_that( + is.character(filename), + length(filename) == 1, + msg = "`filename` must be a length 1 character" + ) + if (nchar(filename) > 1000) { + return(FALSE) + } + if (is.na(filename)) { + return(FALSE) + } + return(file.exists(filename) & !dir.exists(filename)) } @@ -227,16 +225,16 @@ is_file <- function(filename = NULL) { #' file in the package directory or the stan code as a string. #' @export read_stan <- function(string) { - local_inst_file <- file.path("inst", "stan", string) - system_file <- system.file(file.path("stan", string), package = "jmpost") - local_file <- string - files <- c(local_file, local_inst_file, system_file) - for (fi in files) { - if (is_file(fi)) { - return(readLines(fi)) - } + local_inst_file <- file.path("inst", "stan", string) + system_file <- system.file(file.path("stan", string), package = "jmpost") + local_file <- string + files <- c(local_file, local_inst_file, system_file) + for (fi in files) { + if (is_file(fi)) { + return(readLines(fi)) } - return(string) + } + return(string) } @@ -248,20 +246,19 @@ as_stan_file <- function( parameters = "", transformed_parameters = "", model = "", - generated_quantities = "" -) { - block_strings <- lapply( - names(STAN_BLOCKS), - function(id) { - char <- get(id) - if (any(nchar(char) >= 1)) { - return(sprintf("%s {\n%s\n}\n\n", STAN_BLOCKS[[id]], paste(char, collapse = "\n"))) - } else { - return("") - } - } - ) - return(paste0(block_strings, collapse = "")) + generated_quantities = "") { + block_strings <- lapply( + names(STAN_BLOCKS), + function(id) { + char <- get(id) + if (any(nchar(char) >= 1)) { + return(sprintf("%s {\n%s\n}\n\n", STAN_BLOCKS[[id]], paste(char, collapse = "\n"))) + } else { + return("") + } + } + ) + return(paste0(block_strings, collapse = "")) } @@ -283,43 +280,42 @@ as_stan_file <- function( # ``` #' @import stringr as_stan_fragments <- function(x) { - code <- unlist(stringr::str_split(x, "\n")) - results <- list() - target <- NULL - for (line in code) { - for (block in names(STAN_BLOCKS)) { - regex <- sprintf("^%s *\\{ *$", STAN_BLOCKS[[block]]) - if (stringr::str_detect(line, regex)) { - target <- block - line <- NULL - break - } - } - if(!is.null(target)) { - results[[target]] <- c(results[[target]], line) - } + code <- unlist(stringr::str_split(x, "\n")) + results <- list() + target <- NULL + for (line in code) { + for (block in names(STAN_BLOCKS)) { + regex <- sprintf("^%s *\\{ *$", STAN_BLOCKS[[block]]) + if (stringr::str_detect(line, regex)) { + target <- block + line <- NULL + break + } } - - # Remove trailing "}" - for (block in names(results)) { - block_length <- length(results[[block]]) - entry <- block_length - while (entry >= 0) { - line <- results[[block]][[entry]] - if (stringr::str_detect(line, "^ *\\} *$")) { - break - } - entry <- entry - 1 - } - results[[block]] <- results[[block]][-seq(entry, block_length)] + if (!is.null(target)) { + results[[target]] <- c(results[[target]], line) } + } + + # Remove trailing "}" + for (block in names(results)) { + block_length <- length(results[[block]]) + entry <- block_length + while (entry >= 0) { + line <- results[[block]][[entry]] + if (stringr::str_detect(line, "^ *\\} *$")) { + break + } + entry <- entry - 1 + } + results[[block]] <- results[[block]][-seq(entry, block_length)] + } - # Add missings - for (block in names(STAN_BLOCKS)) { - if(is.null(results[[block]])) { - results[[block]] <- "" - } + # Add missings + for (block in names(STAN_BLOCKS)) { + if (is.null(results[[block]])) { + results[[block]] <- "" } - return(results) + } + return(results) } - diff --git a/R/SurvivalExponential.R b/R/SurvivalExponential.R index 6995f820..c9a966c6 100755 --- a/R/SurvivalExponential.R +++ b/R/SurvivalExponential.R @@ -1,26 +1,22 @@ - - #' @include SurvivalModel.R NULL .SurvivalExponential <- setClass( - Class = "SurvivalExponential", - contains = "SurvivalModel" + Class = "SurvivalExponential", + contains = "SurvivalModel" ) #' @export SurvivalExponential <- function( lambda = prior_gamma(2, 5), - beta = prior_normal(0, 5) -) { - .SurvivalExponential( - SurvivalModel( - stan = StanModule("sm-exponential/model.stan"), - parameters = ParameterList( - Parameter(name = "sm_exp_lambda", prior = lambda, size = 1), - Parameter(name = "beta_os_cov", prior = beta, size = "p_os_cov_design") - ) - ) + beta = prior_normal(0, 5)) { + .SurvivalExponential( + SurvivalModel( + stan = StanModule("sm-exponential/model.stan"), + parameters = ParameterList( + Parameter(name = "sm_exp_lambda", prior = lambda, size = 1), + Parameter(name = "beta_os_cov", prior = beta, size = "p_os_cov_design") + ) ) + ) } - diff --git a/R/SurvivalLoglogistic.R b/R/SurvivalLoglogistic.R index 195feea9..297f3474 100644 --- a/R/SurvivalLoglogistic.R +++ b/R/SurvivalLoglogistic.R @@ -1,32 +1,23 @@ - - - #' @include SurvivalModel.R NULL .SurvivalLogLogistic <- setClass( - Class = "SurvivalLogLogistic", - contains = "SurvivalModel" + Class = "SurvivalLogLogistic", + contains = "SurvivalModel" ) #' @export -SurvivalLogLogistic <- function( - lambda = prior_lognormal(0, 5, init = 1 / 200), - p = prior_gamma(2, 5, init = 0.5), - beta = prior_normal(0, 5) -) { - .SurvivalLogLogistic( - SurvivalModel( - stan = StanModule("sm-loglogistic/model.stan"), - parameters = ParameterList( - Parameter(name = "sm_logl_lambda", prior = lambda, size = 1), - Parameter(name = "sm_logl_p", prior = p, size = 1), - Parameter(name = "beta_os_cov", prior = beta, size = "p_os_cov_design") - ) - ) +SurvivalLogLogistic <- function(lambda = prior_lognormal(0, 5, init = 1 / 200), + p = prior_gamma(2, 5, init = 0.5), + beta = prior_normal(0, 5)) { + .SurvivalLogLogistic( + SurvivalModel( + stan = StanModule("sm-loglogistic/model.stan"), + parameters = ParameterList( + Parameter(name = "sm_logl_lambda", prior = lambda, size = 1), + Parameter(name = "sm_logl_p", prior = p, size = 1), + Parameter(name = "beta_os_cov", prior = beta, size = "p_os_cov_design") + ) ) + ) } - - - - diff --git a/R/SurvivalModel.R b/R/SurvivalModel.R index 363b2588..f599bd2a 100755 --- a/R/SurvivalModel.R +++ b/R/SurvivalModel.R @@ -1,26 +1,25 @@ - #' @include StanModel.R NULL .SurvivalModel <- setClass( - Class = "SurvivalModel", - contains = "StanModel" + Class = "SurvivalModel", + contains = "StanModel" ) #' @export SurvivalModel <- function(stan = StanModule(), parameters = ParameterList(), ...) { - base_stan <- paste0(read_stan("base/survival.stan"), collapse = "\n") - stan_full <- jinjar::render( - .x = base_stan, - stan = add_missing_stan_blocks(as.list(stan)) - ) - .SurvivalModel( - StanModel( - stan = StanModule(stan_full), - parameters = parameters, - ... - ) + base_stan <- paste0(read_stan("base/survival.stan"), collapse = "\n") + stan_full <- jinjar::render( + .x = base_stan, + stan = add_missing_stan_blocks(as.list(stan)) + ) + .SurvivalModel( + StanModel( + stan = StanModule(stan_full), + parameters = parameters, + ... ) + ) } diff --git a/R/SurvivalWeibullPH.R b/R/SurvivalWeibullPH.R index f244dd1c..ab7da43b 100755 --- a/R/SurvivalWeibullPH.R +++ b/R/SurvivalWeibullPH.R @@ -1,10 +1,9 @@ - #' @include SurvivalModel.R NULL .SurvivalWeibullPH <- setClass( - Class = "SurvivalWeibullPH", - contains = "SurvivalModel" + Class = "SurvivalWeibullPH", + contains = "SurvivalModel" ) @@ -12,17 +11,15 @@ NULL SurvivalWeibullPH <- function( lambda = prior_gamma(2, 0.5), gamma = prior_gamma(2, 0.5), - beta = prior_normal(0, 5) -) { - .SurvivalWeibullPH( - SurvivalModel( - stan = StanModule(x = "sm-weibull-ph/model.stan"), - parameters = ParameterList( - Parameter(name = "sm_weibull_ph_lambda", prior = lambda, size = 1), - Parameter(name = "sm_weibull_ph_gamma", prior = gamma, size = 1), - Parameter(name = "beta_os_cov", prior = beta, size = "p_os_cov_design") - ) - ) + beta = prior_normal(0, 5)) { + .SurvivalWeibullPH( + SurvivalModel( + stan = StanModule(x = "sm-weibull-ph/model.stan"), + parameters = ParameterList( + Parameter(name = "sm_weibull_ph_lambda", prior = lambda, size = 1), + Parameter(name = "sm_weibull_ph_gamma", prior = gamma, size = 1), + Parameter(name = "beta_os_cov", prior = beta, size = "p_os_cov_design") + ) ) + ) } - diff --git a/R/defaults.R b/R/defaults.R index afab327e..df516f3a 100755 --- a/R/defaults.R +++ b/R/defaults.R @@ -1,4 +1,3 @@ - #' @include generics.R #' @include Link.R #' @include LongitudinalModel.R @@ -15,23 +14,23 @@ NULL # setMethod( - "addLink", - signature = c("LongitudinalModel", "NULL"), - definition = function(x, y, ...) x + "addLink", + signature = c("LongitudinalModel", "NULL"), + definition = function(x, y, ...) x ) setMethod( - "addLink", - signature = c("NULL", "Link"), - definition = function(x, y, ...) { - stop("No Longitudinal model has been defined for the link function to be combined with") - } + "addLink", + signature = c("NULL", "Link"), + definition = function(x, y, ...) { + stop("No Longitudinal model has been defined for the link function to be combined with") + } ) setMethod( - "addLink", - signature = c("NULL", "NULL"), - definition = function(x, y, ...) NULL + "addLink", + signature = c("NULL", "NULL"), + definition = function(x, y, ...) NULL ) @@ -44,9 +43,9 @@ setMethod( setMethod( - f = "getParameters", - signature = "NULL", - definition = function(object) NULL + f = "getParameters", + signature = "NULL", + definition = function(object) NULL ) @@ -59,50 +58,49 @@ setMethod( setMethod( - "merge", - signature = c("StanModel", "NULL"), - definition = function(x, y, ...) x@stan + "merge", + signature = c("StanModel", "NULL"), + definition = function(x, y, ...) x@stan ) setMethod( - "merge", - signature = c("NULL", "StanModel"), - definition = function(x, y, ...) y@stan + "merge", + signature = c("NULL", "StanModel"), + definition = function(x, y, ...) y@stan ) setMethod( - "merge", - signature = c("StanModule", "NULL"), - definition = function(x, y, ...) x + "merge", + signature = c("StanModule", "NULL"), + definition = function(x, y, ...) x ) setMethod( - "merge", - signature = c("NULL", "StanModule"), - definition = function(x, y, ...) y + "merge", + signature = c("NULL", "StanModule"), + definition = function(x, y, ...) y ) setMethod( - "merge", - signature = c("ParameterList", "NULL"), - definition = function(x, y, ...) x + "merge", + signature = c("ParameterList", "NULL"), + definition = function(x, y, ...) x ) setMethod( - "merge", - signature = c("NULL", "ParameterList"), - definition = function(x, y, ...) y + "merge", + signature = c("NULL", "ParameterList"), + definition = function(x, y, ...) y ) setMethod( - "merge", - signature = c("NULL", "NULL"), - definition = function(x, y, ...) NULL + "merge", + signature = c("NULL", "NULL"), + definition = function(x, y, ...) NULL ) - diff --git a/R/generics.R b/R/generics.R index 6cc86663..6bf448d6 100755 --- a/R/generics.R +++ b/R/generics.R @@ -1,4 +1,3 @@ - # "missing" = no argument provided # "NULL" = explicit NULL setClassUnion("empty", c("missing", "NULL")) @@ -11,8 +10,8 @@ setClassUnion("empty", c("missing", "NULL")) #' @param object TODO #' @param ... TODO setGeneric( - name = "merge", - def = function(x, y, ...) standardGeneric("merge") + name = "merge", + def = function(x, y, ...) standardGeneric("merge") ) @@ -24,8 +23,8 @@ setGeneric( #' @param object TODO #' @param ... TODO setGeneric( - name = "addLink", - def = function(x, y, ...) standardGeneric("addLink") + name = "addLink", + def = function(x, y, ...) standardGeneric("addLink") ) @@ -36,8 +35,8 @@ setGeneric( #' @param object TODO #' @param file_path TODO setGeneric( - name = "write_stan", - def = function(x, file_path) standardGeneric("write_stan") + name = "write_stan", + def = function(x, file_path) standardGeneric("write_stan") ) @@ -49,8 +48,8 @@ setGeneric( #' @param exe_file TODO #' @export setGeneric( - name = "compileStanModel", - def = function(object, exe_file) standardGeneric("compileStanModel") + name = "compileStanModel", + def = function(object, exe_file) standardGeneric("compileStanModel") ) @@ -63,8 +62,8 @@ setGeneric( #' @param exe_file TODO #' @export setGeneric( - name = "sampleStanModel", - def = function(object, ..., exe_file) standardGeneric("sampleStanModel") + name = "sampleStanModel", + def = function(object, ..., exe_file) standardGeneric("sampleStanModel") ) @@ -77,8 +76,8 @@ setGeneric( #' @param object TODO #' @param ... TODO setGeneric( - name = "as.StanModule", - def = function(object, ...) standardGeneric("as.StanModule") + name = "as.StanModule", + def = function(object, ...) standardGeneric("as.StanModule") ) @@ -89,8 +88,8 @@ setGeneric( #' @param object TODO #' @param ... TODO setGeneric( - name = "getParameters", - def = function(object, ...) standardGeneric("getParameters") + name = "getParameters", + def = function(object, ...) standardGeneric("getParameters") ) @@ -104,34 +103,34 @@ setGeneric( #' #' @keywords internal setGeneric( - name = "extractVariableNames", - def = function(object, ...) standardGeneric("extractVariableNames") + name = "extractVariableNames", + def = function(object, ...) standardGeneric("extractVariableNames") ) #' `initialValues` -#' +#' #' @param object TODO #' @param ... TODO #' @return a `list` of initial values to be passed to the stan sampler -#' +#' #' @keywords internal setGeneric( - name = "initialValues", - def = function(object, ...) standardGeneric("initialValues") + name = "initialValues", + def = function(object, ...) standardGeneric("initialValues") ) #' `size` -#' +#' #' @param object TODO #' @param ... TODO #' @return a `list` of parameter sizes #' @keywords internal setGeneric( - name = "size", - def = function(object, ...) standardGeneric("size") + name = "size", + def = function(object, ...) standardGeneric("size") ) diff --git a/R/simulations.R b/R/simulations.R index 91ace0ea..f7764458 100755 --- a/R/simulations.R +++ b/R/simulations.R @@ -1,17 +1,15 @@ - - #' Construct a Variance-Covariance Matrix from Standard Deviations and Correlations #' #' This function creates a variance-covariance matrix based on the input standard deviations #' and correlations. Note that the values may be altered using [Matrix::nearPD()] -#' in order to ensure that it is a valid positive semi-definite matrix. +#' in order to ensure that it is a valid positive semi-definite matrix. #' #' @param sd A numeric vector containing the standard deviations of the variables. #' @param cor A numeric vector containing the pairwise correlations between the variables. #' The vector should be in the order of the upper triangular part of the matrix, excluding #' the diagonal elements. e.g. for a 3x3 matrix c(0.2, 0.3, 0.4) would be the values for #' x_12, x_13 and x_23 of the correlation matrix respectively -#' +#' #' #' @return A symmetric matrix representing the variance-covariance matrix calculated from #' the input standard deviations and correlations. The matrix is guaranteed to be positive @@ -24,35 +22,35 @@ #' #' @export as_vcov <- function(sd, cor) { - x <- diag(rep(1, length(sd))) - x[upper.tri(x)] <- cor - x <- t(x) - x[upper.tri(x)] <- cor - res <- diag(sd) %*% x %*% diag(sd) - res <- as.matrix(Matrix::nearPD(res)$mat) - assertthat::assert_that(isSymmetric(res)) - dimnames(res) <- NULL - return(res) + x <- diag(rep(1, length(sd))) + x[upper.tri(x)] <- cor + x <- t(x) + x[upper.tri(x)] <- cor + res <- diag(sd) %*% x %*% diag(sd) + res <- as.matrix(Matrix::nearPD(res)$mat) + assertthat::assert_that(isSymmetric(res)) + dimnames(res) <- NULL + return(res) } #' @export sld <- function(time, b, s, g, phi) { - b * (phi * exp(-s * time) + (1 - phi) * exp(g * time)) + b * (phi * exp(-s * time) + (1 - phi) * exp(g * time)) } ttg <- function(time, b, s, g, phi) { - t1 <- (log(s * phi / (g * (1 - phi))) / (g + s)) - t1[t1 <= 0] <- 0 - return(t1) + t1 <- (log(s * phi / (g * (1 - phi))) / (g + s)) + t1[t1 <= 0] <- 0 + return(t1) } -dsld <- function(time, b , s, g, phi) { - t1 <- (1 - phi) * g * exp(g * time) - t2 <- phi * s * exp(-s * time) - return(b * (t1 - t2)) +dsld <- function(time, b, s, g, phi) { + t1 <- (1 - phi) * g * exp(g * time) + t2 <- phi * s * exp(-s * time) + return(b * (t1 - t2)) } @@ -74,47 +72,45 @@ sim_lm_gsf <- function( omega_phi = 0.75, link_dsld = 0, link_ttg = 0, - .debug = FALSE -) { - function(lm_base) { - - assertthat::assert_that( - length(unique(lm_base$study)) == 1, - length(mu_b) == 1, - length(sigma) == 1, - length(mu_s) == length(unique(lm_base$arm)), - length(mu_s) == length(mu_g), - length(mu_s) == length(mu_phi), - length(c(eta_b_sigma, eta_g_sigma, eta_phi_sigma, eta_s_sigma)) == 4, - length(c(omega_b, omega_s, omega_g, omega_phi)) == 4 - ) - - baseline_covs <- lm_base |> - dplyr::distinct(pt, arm, study) |> - dplyr::mutate(arm_n = as.numeric(factor(as.character(arm)))) |> - dplyr::mutate(eta_b = rnorm(dplyr::n(), 0, eta_b_sigma)) |> - dplyr::mutate(eta_s = rnorm(dplyr::n(), 0, eta_s_sigma)) |> - dplyr::mutate(eta_g = rnorm(dplyr::n(), 0, eta_g_sigma)) |> - dplyr::mutate(eta_phi = rnorm(dplyr::n(), 0, eta_phi_sigma)) |> - dplyr::mutate(psi_b = exp(log(mu_b) + eta_b * omega_b)) |> - dplyr::mutate(psi_s = exp(log(mu_s[arm_n]) + eta_s * omega_s)) |> - dplyr::mutate(psi_g = exp(log(mu_g[arm_n]) + eta_g * omega_g)) |> - dplyr::mutate(psi_phi = plogis(qlogis(mu_phi[arm_n]) + eta_phi * omega_phi)) - - lm_dat <- lm_base |> - dplyr::select(-study, -arm) |> - dplyr::left_join(baseline_covs, by = "pt") |> - dplyr::mutate(mu_sld = sld(time, psi_b, psi_s, psi_g, psi_phi)) |> - dplyr::mutate(dsld = dsld(time, psi_b, psi_s, psi_g, psi_phi)) |> - dplyr::mutate(ttg = ttg(time, psi_b, psi_s, psi_g, psi_phi)) |> - dplyr::mutate(sld = rnorm(dplyr::n(), mu_sld, mu_sld * sigma)) |> - dplyr::mutate(log_haz_link = link_dsld * dsld + link_ttg * ttg) - - if (!.debug) { - lm_dat <- lm_dat |> dplyr::select(pt, time, sld, log_haz_link, study, arm) - } - return(lm_dat) + .debug = FALSE) { + function(lm_base) { + assertthat::assert_that( + length(unique(lm_base$study)) == 1, + length(mu_b) == 1, + length(sigma) == 1, + length(mu_s) == length(unique(lm_base$arm)), + length(mu_s) == length(mu_g), + length(mu_s) == length(mu_phi), + length(c(eta_b_sigma, eta_g_sigma, eta_phi_sigma, eta_s_sigma)) == 4, + length(c(omega_b, omega_s, omega_g, omega_phi)) == 4 + ) + + baseline_covs <- lm_base |> + dplyr::distinct(pt, arm, study) |> + dplyr::mutate(arm_n = as.numeric(factor(as.character(arm)))) |> + dplyr::mutate(eta_b = rnorm(dplyr::n(), 0, eta_b_sigma)) |> + dplyr::mutate(eta_s = rnorm(dplyr::n(), 0, eta_s_sigma)) |> + dplyr::mutate(eta_g = rnorm(dplyr::n(), 0, eta_g_sigma)) |> + dplyr::mutate(eta_phi = rnorm(dplyr::n(), 0, eta_phi_sigma)) |> + dplyr::mutate(psi_b = exp(log(mu_b) + eta_b * omega_b)) |> + dplyr::mutate(psi_s = exp(log(mu_s[arm_n]) + eta_s * omega_s)) |> + dplyr::mutate(psi_g = exp(log(mu_g[arm_n]) + eta_g * omega_g)) |> + dplyr::mutate(psi_phi = plogis(qlogis(mu_phi[arm_n]) + eta_phi * omega_phi)) + + lm_dat <- lm_base |> + dplyr::select(-study, -arm) |> + dplyr::left_join(baseline_covs, by = "pt") |> + dplyr::mutate(mu_sld = sld(time, psi_b, psi_s, psi_g, psi_phi)) |> + dplyr::mutate(dsld = dsld(time, psi_b, psi_s, psi_g, psi_phi)) |> + dplyr::mutate(ttg = ttg(time, psi_b, psi_s, psi_g, psi_phi)) |> + dplyr::mutate(sld = rnorm(dplyr::n(), mu_sld, mu_sld * sigma)) |> + dplyr::mutate(log_haz_link = link_dsld * dsld + link_ttg * ttg) + + if (!.debug) { + lm_dat <- lm_dat |> dplyr::select(pt, time, sld, log_haz_link, study, arm) } + return(lm_dat) + } } @@ -125,189 +121,185 @@ sim_lm_random_slope <- function( slope_sigma = 0.5, sigma = 2, phi = 0.1, - .debug = FALSE -) { - function(lm_base) { - - assertthat::assert_that( - length(slope_mu) == 1 | length(slope_mu) == length(unique(lm_base$arm)), - msg = "slope_mu should either be length 1 or equal to the length of n_arm" - ) - - if (length(slope_mu) == 1) { - slope_mu <- rep(slope_mu, length(unique(lb_base$arm))) - } - - rs_baseline <- lm_base |> - dplyr::distinct(pt, arm) |> - dplyr::mutate(slope_ind = rnorm(dplyr::n(), slope_mu[as.numeric(arm)], sd = slope_sigma)) |> - dplyr::select(-arm) - - lm_dat <- lm_base |> - dplyr::mutate(err = rnorm(dplyr::n(), 0, sigma)) |> - dplyr::left_join(rs_baseline, by = "pt") |> - dplyr::mutate(sld = intercept + slope_ind * time + err) |> - dplyr::mutate(log_haz_link = slope_ind * phi) - - if ( ! .debug) { - lm_dat <- lm_dat |> dplyr::select(pt, time, sld, log_haz_link, study, arm) - } - return(lm_dat) + .debug = FALSE) { + function(lm_base) { + assertthat::assert_that( + length(slope_mu) == 1 | length(slope_mu) == length(unique(lm_base$arm)), + msg = "slope_mu should either be length 1 or equal to the length of n_arm" + ) + + if (length(slope_mu) == 1) { + slope_mu <- rep(slope_mu, length(unique(lb_base$arm))) } + + rs_baseline <- lm_base |> + dplyr::distinct(pt, arm) |> + dplyr::mutate(slope_ind = rnorm(dplyr::n(), slope_mu[as.numeric(arm)], sd = slope_sigma)) |> + dplyr::select(-arm) + + lm_dat <- lm_base |> + dplyr::mutate(err = rnorm(dplyr::n(), 0, sigma)) |> + dplyr::left_join(rs_baseline, by = "pt") |> + dplyr::mutate(sld = intercept + slope_ind * time + err) |> + dplyr::mutate(log_haz_link = slope_ind * phi) + + if (!.debug) { + lm_dat <- lm_dat |> dplyr::select(pt, time, sld, log_haz_link, study, arm) + } + return(lm_dat) + } } #' @export sim_os_weibull <- function(lambda, gamma) { - function(time) { - log(lambda) + log(gamma) + (gamma - 1) * log(time) - } + function(time) { + log(lambda) + log(gamma) + (gamma - 1) * log(time) + } } #' @export sim_os_exponential <- function(lambda) { - function(time) { - log(lambda) - } + function(time) { + log(lambda) + } } #' @export -sim_os_loglogistic <- function(lambda, p){ - function(time) { - c1 <- log(lambda) + log(p) + (p - 1) * (log(lambda) + log(time)) - c2 <- log(1 + (lambda * time)^p) - return(c1 - c2) - } +sim_os_loglogistic <- function(lambda, p) { + function(time) { + c1 <- log(lambda) + log(p) + (p - 1) * (log(lambda) + log(time)) + c2 <- log(1 + (lambda * time)^p) + return(c1 - c2) + } } get_timepoints <- function(x) { - assertthat::assert_that(length(x) == length(unique(x))) - x_ord <- x[order(x)] - x_no_neg <- x_ord[x_ord > 0] - - bound_lower <- c(0, x_no_neg[-length(x_no_neg)]) - bound_upper <- x_no_neg - bound_width <- bound_upper - bound_lower - eval_points <- bound_upper - tibble::tibble( - lower = bound_lower, - upper = bound_upper, - time = eval_points, - width = bound_width - ) + assertthat::assert_that(length(x) == length(unique(x))) + x_ord <- x[order(x)] + x_no_neg <- x_ord[x_ord > 0] + + bound_lower <- c(0, x_no_neg[-length(x_no_neg)]) + bound_upper <- x_no_neg + bound_width <- bound_upper - bound_lower + eval_points <- bound_upper + tibble::tibble( + lower = bound_lower, + upper = bound_upper, + time = eval_points, + width = bound_width + ) } #' @export simulate_joint_data <- function( - n_arm = c(50, 80), # Number of arms and number of subjects per arm + n_arm = c(50, 80), # Number of arms and number of subjects per arm times = 1:2000, lambda_cen = 1 / 3, beta_cont = 0.2, beta_cat = c("A" = 0, "B" = -0.4, "C" = 0.2), lm_fun, - os_fun -) { - - n <- sum(n_arm) - u_pts <- sprintf("pt_%05i", seq_len(n)) - bounds <- get_timepoints(times) - - ARMS <- paste0("Group-", seq_along(n_arm)) - - os_baseline <- dplyr::tibble(pt = u_pts) |> - dplyr::mutate(cov_cont = rnorm(n)) |> - dplyr::mutate(cov_cat = factor( - sample(c("A", "B", "C"), replace = TRUE, size = n), - levels = c("A", "B", "C") - )) |> - dplyr::mutate(log_haz_cov = cov_cont * beta_cont + beta_cat[cov_cat]) |> - dplyr::mutate(survival = runif(n)) |> - dplyr::mutate(chazard_limit = -log(survival)) |> - dplyr::mutate(time_cen = rexp(n, lambda_cen)) |> - dplyr::mutate(study = factor("Study-1")) |> - dplyr::mutate(arm = factor(rep(ARMS, times = n_arm), levels = ARMS)) - - time_dat <- expand.grid( - time = as.double(bounds$time), - pt = u_pts, - stringsAsFactors = FALSE - ) |> - tibble::as_tibble() |> - dplyr::mutate(width = rep(bounds$width, times = n)) - - time_dat_baseline <- time_dat |> - dplyr::left_join(os_baseline, by = "pt") - - lm_dat <- time_dat_baseline |> - dplyr::select(pt, time, arm, study) |> - lm_fun() - - assert_that( - all(lm_dat$pt == time_dat_baseline$pt), - all(lm_dat$time == time_dat_baseline$time), - msg = "The longitudinal dataset must be sorted by pt, time" - ) - - os_dat_chaz <- time_dat_baseline |> - dplyr::mutate(log_haz_link = lm_dat$log_haz_link) |> # only works if lm_dat is sorted pt, time - dplyr::mutate(log_bl_haz = os_fun(time)) |> - dplyr::mutate(hazard_instant = exp(log_bl_haz + log_haz_cov + log_haz_link)) |> - dplyr::mutate(hazard_interval = hazard_instant * width) |> - dplyr::group_by(pt) |> - dplyr::mutate(chazard = cumsum(hazard_interval)) |> - dplyr::ungroup() - - os_dat <- os_dat_chaz|> - dplyr::filter(chazard_limit <= chazard) |> - dplyr::select(pt, time, cov_cont, cov_cat, time_cen, study, arm) |> - dplyr::group_by(pt) |> - dplyr::slice(1) |> - dplyr::ungroup() |> - dplyr::mutate(event = dplyr::if_else(time <= time_cen, 1, 0)) |> - dplyr::mutate(time = dplyr::if_else(event == 1, time, time_cen)) |> - dplyr::select(pt, time, cov_cont, cov_cat, event, study, arm) - - os_dat_censor <- os_dat_chaz |> - dplyr::group_by(pt) |> - dplyr::slice(1) |> - dplyr::mutate(time = max(times)) |> - dplyr::mutate(event = 0) |> - dplyr::select(pt, time, cov_cont, cov_cat, event, study, arm) |> - dplyr::filter(!pt %in% os_dat$pt) - - if (!nrow(os_dat_censor)== 0 ) { - message(sprintf("%i patients did not die before max(times)", nrow(os_dat_censor))) - } - - os_dat_complete <- os_dat |> - dplyr::bind_rows(os_dat_censor) |> - dplyr::arrange(pt) - - os_time <- rep(os_dat_complete$time, each = length(bounds$time)) - - assert_that( - length(os_time) == nrow(lm_dat), - all(lm_dat$pt == rep(os_dat_complete$pt, each = length(bounds$time))) - ) - - lm_dat2 <- lm_dat |> - dplyr::mutate(os_time = os_time) |> - dplyr::mutate(observed = (time <= os_time)) |> - dplyr::select(-os_time, -log_haz_link) - - assertthat::assert_that( - length(unique(os_dat_complete$pt)) == n, - n == nrow(os_dat_complete), - all(os_dat_complete$time >= 0), - all(os_dat_complete$event %in% c(0, 1)), - nrow(lm_dat2) == n * length(bounds$time), - all(!is.na(lm_dat2$observed)) - ) - return(list(os = os_dat_complete, lm = lm_dat2)) + os_fun) { + n <- sum(n_arm) + u_pts <- sprintf("pt_%05i", seq_len(n)) + bounds <- get_timepoints(times) + + ARMS <- paste0("Group-", seq_along(n_arm)) + + os_baseline <- dplyr::tibble(pt = u_pts) |> + dplyr::mutate(cov_cont = rnorm(n)) |> + dplyr::mutate(cov_cat = factor( + sample(c("A", "B", "C"), replace = TRUE, size = n), + levels = c("A", "B", "C") + )) |> + dplyr::mutate(log_haz_cov = cov_cont * beta_cont + beta_cat[cov_cat]) |> + dplyr::mutate(survival = runif(n)) |> + dplyr::mutate(chazard_limit = -log(survival)) |> + dplyr::mutate(time_cen = rexp(n, lambda_cen)) |> + dplyr::mutate(study = factor("Study-1")) |> + dplyr::mutate(arm = factor(rep(ARMS, times = n_arm), levels = ARMS)) + + time_dat <- expand.grid( + time = as.double(bounds$time), + pt = u_pts, + stringsAsFactors = FALSE + ) |> + tibble::as_tibble() |> + dplyr::mutate(width = rep(bounds$width, times = n)) + + time_dat_baseline <- time_dat |> + dplyr::left_join(os_baseline, by = "pt") + + lm_dat <- time_dat_baseline |> + dplyr::select(pt, time, arm, study) |> + lm_fun() + + assert_that( + all(lm_dat$pt == time_dat_baseline$pt), + all(lm_dat$time == time_dat_baseline$time), + msg = "The longitudinal dataset must be sorted by pt, time" + ) + + os_dat_chaz <- time_dat_baseline |> + dplyr::mutate(log_haz_link = lm_dat$log_haz_link) |> # only works if lm_dat is sorted pt, time + dplyr::mutate(log_bl_haz = os_fun(time)) |> + dplyr::mutate(hazard_instant = exp(log_bl_haz + log_haz_cov + log_haz_link)) |> + dplyr::mutate(hazard_interval = hazard_instant * width) |> + dplyr::group_by(pt) |> + dplyr::mutate(chazard = cumsum(hazard_interval)) |> + dplyr::ungroup() + + os_dat <- os_dat_chaz |> + dplyr::filter(chazard_limit <= chazard) |> + dplyr::select(pt, time, cov_cont, cov_cat, time_cen, study, arm) |> + dplyr::group_by(pt) |> + dplyr::slice(1) |> + dplyr::ungroup() |> + dplyr::mutate(event = dplyr::if_else(time <= time_cen, 1, 0)) |> + dplyr::mutate(time = dplyr::if_else(event == 1, time, time_cen)) |> + dplyr::select(pt, time, cov_cont, cov_cat, event, study, arm) + + os_dat_censor <- os_dat_chaz |> + dplyr::group_by(pt) |> + dplyr::slice(1) |> + dplyr::mutate(time = max(times)) |> + dplyr::mutate(event = 0) |> + dplyr::select(pt, time, cov_cont, cov_cat, event, study, arm) |> + dplyr::filter(!pt %in% os_dat$pt) + + if (!nrow(os_dat_censor) == 0) { + message(sprintf("%i patients did not die before max(times)", nrow(os_dat_censor))) + } + + os_dat_complete <- os_dat |> + dplyr::bind_rows(os_dat_censor) |> + dplyr::arrange(pt) + + os_time <- rep(os_dat_complete$time, each = length(bounds$time)) + + assert_that( + length(os_time) == nrow(lm_dat), + all(lm_dat$pt == rep(os_dat_complete$pt, each = length(bounds$time))) + ) + + lm_dat2 <- lm_dat |> + dplyr::mutate(os_time = os_time) |> + dplyr::mutate(observed = (time <= os_time)) |> + dplyr::select(-os_time, -log_haz_link) + + assertthat::assert_that( + length(unique(os_dat_complete$pt)) == n, + n == nrow(os_dat_complete), + all(os_dat_complete$time >= 0), + all(os_dat_complete$event %in% c(0, 1)), + nrow(lm_dat2) == n * length(bounds$time), + all(!is.na(lm_dat2$observed)) + ) + return(list(os = os_dat_complete, lm = lm_dat2)) } diff --git a/R/utilities.R b/R/utilities.R index cec8e3e9..8e0dfd2b 100644 --- a/R/utilities.R +++ b/R/utilities.R @@ -1,54 +1,51 @@ - - - #' `get_missing_rownumbers` -#' +#' #' Returns which row numbers contain at least 1 missing observation any variables -#' +#' #' @param df A data.frame #' @param formula A formula specifying which variables to inspect for missingness, if null #' all variables are considered -#' +#' #' @keywords internal get_missing_rownumbers <- function(df, formula = NULL) { - if (is.null(formula)) { - formula <- ~ . - } - mdf <- model.frame(formula, data = df, na.action = na.pass) - missing_rows <- which(!complete.cases(mdf)) - return(missing_rows) + if (is.null(formula)) { + formula <- ~. + } + mdf <- model.frame(formula, data = df, na.action = na.pass) + missing_rows <- which(!complete.cases(mdf)) + return(missing_rows) } #' `remove_missing_rows` -#' +#' #' Removes any rows from a dataset that contain missing values. Allows users to specify which variables #' to inspect for missing values based on either a formula or a character vector of variable names. -#' +#' #' @param data A data.frame #' @param formula A formula specifying which variables to inspect for missingness #' @param extra_vars A character vector specifying which variables to inspect for missingness -#' +#' #' @returns `data` after removing observations that contain missing values in the required variables -#' Note that additional variables not listed in `formula` or `extra_vars` are not dropped and may +#' Note that additional variables not listed in `formula` or `extra_vars` are not dropped and may #' still contain missing values -#' +#' #' @keywords internal remove_missing_rows <- function(data, formula, extra_vars) { - extra_vars <- paste(extra_vars, collapse = " + ") - formula_update_string <- paste0(". ~ . + ", extra_vars) - formula_update <- update(formula, formula_update_string) - - missing_rows <- get_missing_rownumbers(data, formula_update) - - if (length(missing_rows) == 0 ) { - return(data) - } - message(sprintf( - "Note that %d observations were removed as one of more required variables contained missing values", - length(missing_rows) - )) - return(data[-missing_rows, ]) + extra_vars <- paste(extra_vars, collapse = " + ") + formula_update_string <- paste0(". ~ . + ", extra_vars) + formula_update <- update(formula, formula_update_string) + + missing_rows <- get_missing_rownumbers(data, formula_update) + + if (length(missing_rows) == 0) { + return(data) + } + message(sprintf( + "Note that %d observations were removed as one of more required variables contained missing values", + length(missing_rows) + )) + return(data[-missing_rows, ]) } @@ -75,37 +72,37 @@ remove_missing_rows <- function(data, formula, extra_vars) { #' #' @importFrom assertthat assert_that #' @keywords internal -expand_initial_values <- function(initial_values, sizes){ - assert_that( - is.list(initial_values), - msg = "`initial_values` must be a list" - ) - assert_that( - is.list(sizes), - msg = "`sizes` must be a list" - ) +expand_initial_values <- function(initial_values, sizes) { + assert_that( + is.list(initial_values), + msg = "`initial_values` must be a list" + ) + assert_that( + is.list(sizes), + msg = "`sizes` must be a list" + ) + assert_that( + all(names(sizes) %in% names(initial_values)), + all(names(initial_values) %in% names(sizes)), + msg = "`initial_values` and `sizes` must have identical names" + ) + + # Check for single values in initial_values and replicate them according to sizes + for (name in names(initial_values)) { + if (length(initial_values[[name]]) == 1) { + initial_values[[name]] <- rep(initial_values[[name]], sizes[[name]]) + } + } + + # Check that each element of initial_values has the same length as specified in sizes + for (name in names(initial_values)) { assert_that( - all(names(sizes) %in% names(initial_values)), - all(names(initial_values) %in% names(sizes)), - msg = "`initial_values` and `sizes` must have identical names" + length(initial_values[[name]]) == sizes[[name]], + msg = "length of element in `initial_values` does not match specified size" ) - - # Check for single values in initial_values and replicate them according to sizes - for (name in names(initial_values)) { - if (length(initial_values[[name]]) == 1) { - initial_values[[name]] <- rep(initial_values[[name]], sizes[[name]]) - } - } - - # Check that each element of initial_values has the same length as specified in sizes - for (name in names(initial_values)) { - assert_that( - length(initial_values[[name]]) == sizes[[name]], - msg = "length of element in `initial_values` does not match specified size" - ) - } - - return(initial_values) + } + + return(initial_values) } @@ -132,37 +129,34 @@ expand_initial_values <- function(initial_values, sizes){ #' #' @keywords internal replace_with_lookup <- function(sizes, data) { - - assert_that( is.list(sizes), msg = "`sizes` must be a list") - assert_that( is.list(data), msg = "`data` must be a list") - - for ( idx in seq_along(sizes)) { - val <- sizes[[idx]] - if (is.character(val)) { - assert_that( - length(val) == 1, - msg = "character elements of `sizes` must be length 1" - ) - assert_that( - val %in% names(data), - msg = sprintf("`%s` is not available in `data`") - ) - new_val <- data[[val]] - assert_that( - length(new_val) == 1, - is.numeric(new_val), - msg = "Selected values from data must be length 1 numerics" - ) - sizes[[idx]] <- new_val - } - - assert_that( - is.numeric(sizes[[idx]]), - length(sizes[[idx]]) == 1, - msg = "All elements of `sizes` must be length 1 numerics after lookup" - ) + assert_that(is.list(sizes), msg = "`sizes` must be a list") + assert_that(is.list(data), msg = "`data` must be a list") + + for (idx in seq_along(sizes)) { + val <- sizes[[idx]] + if (is.character(val)) { + assert_that( + length(val) == 1, + msg = "character elements of `sizes` must be length 1" + ) + assert_that( + val %in% names(data), + msg = sprintf("`%s` is not available in `data`") + ) + new_val <- data[[val]] + assert_that( + length(new_val) == 1, + is.numeric(new_val), + msg = "Selected values from data must be length 1 numerics" + ) + sizes[[idx]] <- new_val } - return(sizes) -} - + assert_that( + is.numeric(sizes[[idx]]), + length(sizes[[idx]]) == 1, + msg = "All elements of `sizes` must be length 1 numerics after lookup" + ) + } + return(sizes) +} diff --git a/misc/example_of_use.R b/misc/example_of_use.R index a27185e3..56393c89 100755 --- a/misc/example_of_use.R +++ b/misc/example_of_use.R @@ -19,9 +19,9 @@ devtools::load_all(export_all = FALSE) #### Example 1 - Fully specified model - using the defaults for everything jm <- JointModel( - longitudinal = LongitudinalRandomSlope(), - link = LinkRandomSlope(), - survival = SurvivalWeibullPH() + longitudinal = LongitudinalRandomSlope(), + link = LinkRandomSlope(), + survival = SurvivalWeibullPH() ) @@ -30,21 +30,21 @@ jm <- JointModel( ### Example 2 - Manually specify priors - Fit models independently (no link) jm <- JointModel( - link = LinkNone(), - longitudinal = LongitudinalRandomSlope( - intercept = prior_normal(40, 5), # Just prior - slope_mu = prior_normal(10, 2, init = 30) # Prior and init - ), - survival = SurvivalWeibullPH( - lambda = prior_gamma(0.2, 0.5) - ) + link = LinkNone(), + longitudinal = LongitudinalRandomSlope( + intercept = prior_normal(40, 5), # Just prior + slope_mu = prior_normal(10, 2, init = 30) # Prior and init + ), + survival = SurvivalWeibullPH( + lambda = prior_gamma(0.2, 0.5) + ) ) ### Example 3 - Specify survival model only jm <- JointModel( - survival = SurvivalWeibullPH() + survival = SurvivalWeibullPH() ) @@ -52,7 +52,7 @@ jm <- JointModel( ### Example 4 - Specify longitudinal model only jm <- JointModel( - longitudinal = LongitudinalRandomSlope() + longitudinal = LongitudinalRandomSlope() ) @@ -67,9 +67,9 @@ jm <- JointModel( jm <- JointModel( - longitudinal = LongitudinalRandomSlope(), - link = LinkRandomSlope(), - survival = SurvivalWeibullPH() + longitudinal = LongitudinalRandomSlope(), + link = LinkRandomSlope(), + survival = SurvivalWeibullPH() ) @@ -79,77 +79,77 @@ write_stan(jm, "local/debug.stan") ## Generate Test data with known parameters jlist <- simulate_joint_data( - n_arm = c(500, 500), - times = 1:2000, - lambda_cen = 1 / 9000, - beta_cat = c( - "A" = 0, - "B" = -0.1, - "C" = 0.5 - ), - beta_cont = 0.3, - lm_fun = sim_lm_random_slope( - intercept = 30, - slope_mu = c(1, 2), - slope_sigma = 0.2, - sigma = 3, - phi = 0.1 - ), - os_fun = sim_os_weibull( - lambda = 0.00333, # 1/300 - gamma = 0.97 - ) + n_arm = c(500, 500), + times = 1:2000, + lambda_cen = 1 / 9000, + beta_cat = c( + "A" = 0, + "B" = -0.1, + "C" = 0.5 + ), + beta_cont = 0.3, + lm_fun = sim_lm_random_slope( + intercept = 30, + slope_mu = c(1, 2), + slope_sigma = 0.2, + sigma = 3, + phi = 0.1 + ), + os_fun = sim_os_weibull( + lambda = 0.00333, # 1/300 + gamma = 0.97 + ) ) ## Extract data to individual datasets, reduce longitudinal data to specific time points dat_os <- jlist$os dat_lm <- jlist$lm |> - dplyr::filter(time %in% c(1, 50, 100, 150, 200, 250, 300)) |> - dplyr::arrange(time, pt) + dplyr::filter(time %in% c(1, 50, 100, 150, 200, 250, 300)) |> + dplyr::arrange(time, pt) ## Specify required variables to fit the model to within our dataset jdat <- DataJoint( - survival = DataSurvival( - data = dat_os, - formula = Surv(time, event) ~ cov_cat + cov_cont, - subject = "pt", - arm = "arm", - study = "study" - ), - longitudinal = DataLongitudinal( - data = dat_lm, - formula = sld ~ time, - subject = "pt", - threshold = 5 - ) + survival = DataSurvival( + data = dat_os, + formula = Surv(time, event) ~ cov_cat + cov_cont, + subject = "pt", + arm = "arm", + study = "study" + ), + longitudinal = DataLongitudinal( + data = dat_lm, + formula = sld ~ time, + subject = "pt", + threshold = 5 + ) ) ## Sample from JointModel mp <- sampleStanModel( - jm, - data = jdat, - iter_sampling = 500, - iter_warmup = 500, - chains = 1, - parallel_chains = 1, - exe_file = file.path("local", "full") + jm, + data = jdat, + iter_sampling = 500, + iter_warmup = 500, + chains = 1, + parallel_chains = 1, + exe_file = file.path("local", "full") ) ## Inspect parameters from the model vars <- c( - "lm_rs_intercept", - "lm_rs_slope_mu", - "lm_rs_slope_sigma", - "lm_rs_sigma", - "link_lm_phi", - "sm_weibull_ph_lambda", - "sm_weibull_ph_gamma", - "beta_os_cov" + "lm_rs_intercept", + "lm_rs_slope_mu", + "lm_rs_slope_sigma", + "lm_rs_sigma", + "link_lm_phi", + "sm_weibull_ph_lambda", + "sm_weibull_ph_gamma", + "beta_os_cov" ) diff --git a/misc/tests/gsf-no-link.R b/misc/tests/gsf-no-link.R index 78faebf6..cf906ce2 100644 --- a/misc/tests/gsf-no-link.R +++ b/misc/tests/gsf-no-link.R @@ -18,34 +18,34 @@ devtools::load_all(export_all = FALSE) ## Generate Test data with known parameters jlist <- simulate_joint_data( - n_arm = c(80, 80), - times = seq(0, 4, by = (1/365)/2), - lambda_cen = 1 / 9000, - beta_cat = c( - "A" = 0, - "B" = -0.1, - "C" = 0.5 - ), - beta_cont = 0.3, - lm_fun = sim_lm_gsf( - sigma = 0.003, - mu_s = c(0.2, 0.25), - mu_g = c(0.15, 0.2), - mu_phi = c(0.4, 0.6), - mu_b = 60, - eta_b_sigma = 0.5, - eta_s_sigma = 0.5, - eta_g_sigma = 0.5, - eta_phi_sigma = 0.5, - omega_b = 0.1, - omega_s = 0.1, - omega_g = 0.1, - omega_phi = 0.2 - ), - os_fun = sim_os_weibull( - lambda = 365 * (1/400), - gamma = 1 - ) + n_arm = c(80, 80), + times = seq(0, 4, by = (1 / 365) / 2), + lambda_cen = 1 / 9000, + beta_cat = c( + "A" = 0, + "B" = -0.1, + "C" = 0.5 + ), + beta_cont = 0.3, + lm_fun = sim_lm_gsf( + sigma = 0.003, + mu_s = c(0.2, 0.25), + mu_g = c(0.15, 0.2), + mu_phi = c(0.4, 0.6), + mu_b = 60, + eta_b_sigma = 0.5, + eta_s_sigma = 0.5, + eta_g_sigma = 0.5, + eta_phi_sigma = 0.5, + omega_b = 0.1, + omega_s = 0.1, + omega_g = 0.1, + omega_phi = 0.2 + ), + os_fun = sim_os_weibull( + lambda = 365 * (1 / 400), + gamma = 1 + ) ) @@ -56,9 +56,9 @@ select_times <- c(1, 100, 150, 200, 300, 400, 500, 600, 800, 900) * (1 / 365) # select_times <- seq(1, 2000, by = 30) dat_lm <- jlist$lm |> - dplyr::filter(time %in% select_times) |> - dplyr::arrange(time, pt) |> - dplyr::mutate(time = time) + dplyr::filter(time %in% select_times) |> + dplyr::arrange(time, pt) |> + dplyr::mutate(time = time) # mean(dat_os$time) # mean(dat_os$event) @@ -67,32 +67,28 @@ dat_lm <- jlist$lm |> pnam <- unique(dat_os$pt) |> sample(size = 10) ggplot(data = dat_lm |> dplyr::filter(pt %in% pnam)) + - geom_point(aes(x = time, y = sld, col =pt, group =pt)) + - geom_line(aes(x = time, y = sld, col =pt, group =pt)) + - theme_bw() + geom_point(aes(x = time, y = sld, col = pt, group = pt)) + + geom_line(aes(x = time, y = sld, col = pt, group = pt)) + + theme_bw() jm <- JointModel( - longitudinal = LongitudinalGSF( - - mu_bsld = prior_lognormal(log(60), 2, init = 60), - mu_ks = prior_lognormal(log(0.2), 0.1, init = 0.2), - mu_kg = prior_lognormal(log(0.2), 0.1, init = 0.2), - mu_phi = prior_beta(7, 10, init = 0.4), - - omega_bsld = prior_lognormal(log(0.1), 0.5, init = 0.1), - omega_ks = prior_lognormal(log(0.1), 0.5, init = 0.1), - omega_kg = prior_lognormal(log(0.1), 0.5, init = 0.1), - omega_phi = prior_lognormal(log(0.1), 0.5, init = 0.1), - - sigma = prior_lognormal(log(0.03), 0.1, init = 0.03), - - tilde_bsld = prior_normal(0, 1, init = 0.1), - tilde_ks = prior_normal(0, 1, init = 0.1), - tilde_kg = prior_normal(0, 1, init = 0.1), - tilde_phi = prior_normal(0, 1, init = 0.1) - ) + longitudinal = LongitudinalGSF( + mu_bsld = prior_lognormal(log(60), 2, init = 60), + mu_ks = prior_lognormal(log(0.2), 0.1, init = 0.2), + mu_kg = prior_lognormal(log(0.2), 0.1, init = 0.2), + mu_phi = prior_beta(7, 10, init = 0.4), + omega_bsld = prior_lognormal(log(0.1), 0.5, init = 0.1), + omega_ks = prior_lognormal(log(0.1), 0.5, init = 0.1), + omega_kg = prior_lognormal(log(0.1), 0.5, init = 0.1), + omega_phi = prior_lognormal(log(0.1), 0.5, init = 0.1), + sigma = prior_lognormal(log(0.03), 0.1, init = 0.03), + tilde_bsld = prior_normal(0, 1, init = 0.1), + tilde_ks = prior_normal(0, 1, init = 0.1), + tilde_kg = prior_normal(0, 1, init = 0.1), + tilde_phi = prior_normal(0, 1, init = 0.1) + ) ) @@ -106,45 +102,45 @@ write_stan(jm, "local/debug.stan") ## Prepare data for sampling jdat <- DataJoint( - survival = DataSurvival( - data = dat_os, - formula = Surv(time, event) ~ cov_cat + cov_cont, - subject = "pt", - arm = "arm", - study = "study" - ), - longitudinal = DataLongitudinal( - data = dat_lm, - formula = sld ~ time, - subject = "pt", - threshold = 5 - ) + survival = DataSurvival( + data = dat_os, + formula = Surv(time, event) ~ cov_cat + cov_cont, + subject = "pt", + arm = "arm", + study = "study" + ), + longitudinal = DataLongitudinal( + data = dat_lm, + formula = sld ~ time, + subject = "pt", + threshold = 5 + ) ) ## Sample from JointModel mp <- sampleStanModel( - jm, - data = jdat, - iter_sampling = 500, - iter_warmup = 1000, - chains = 1, - parallel_chains = 1, - exe_file = file.path("local", "full") + jm, + data = jdat, + iter_sampling = 500, + iter_warmup = 1000, + chains = 1, + parallel_chains = 1, + exe_file = file.path("local", "full") ) vars <- c( - "lm_gsf_mu_bsld", # 60 - "lm_gsf_mu_phi", # 0.4 0.6 - "lm_gsf_mu_kg", # 0.15 0.2 - "lm_gsf_mu_ks", # 0.2 0.25 - "lm_gsf_sigma", # 0.003 - "lm_gsf_omega_bsld", # 0.1 - "lm_gsf_omega_kg", # 0.1 - "lm_gsf_omega_phi", # 0.1 - "lm_gsf_omega_ks" # 0.1 + "lm_gsf_mu_bsld", # 60 + "lm_gsf_mu_phi", # 0.4 0.6 + "lm_gsf_mu_kg", # 0.15 0.2 + "lm_gsf_mu_ks", # 0.2 0.25 + "lm_gsf_sigma", # 0.003 + "lm_gsf_omega_bsld", # 0.1 + "lm_gsf_omega_kg", # 0.1 + "lm_gsf_omega_phi", # 0.1 + "lm_gsf_omega_ks" # 0.1 ) @@ -170,14 +166,14 @@ mcmc_pairs(mp$draws(), vars) ### Extract parameters and calculate confidence intervals draws_means <- mp$draws(format = "df") |> - gather() |> - group_by(key) |> - summarise( - mean = mean(value), - sd = sd(value), - lci = quantile(value, 0.025), - uci = quantile(value, 0.975) - ) + gather() |> + group_by(key) |> + summarise( + mean = mean(value), + sd = sd(value), + lci = quantile(value, 0.025), + uci = quantile(value, 0.975) + ) #### Get a list of all named parameters @@ -187,9 +183,9 @@ draws_means |> filter(key == "log_lik[42]") draws_means |> - filter(key %in% vars) |> - mutate(key = factor(key, levels = vars)) |> - arrange(key) + filter(key %in% vars) |> + mutate(key = factor(key, levels = vars)) |> + arrange(key) @@ -199,14 +195,14 @@ draws_means |> devtools::document() devtools::load_all() jm <- JointModel( - survival_model = SurvivalWeibullPH() + survival_model = SurvivalWeibullPH() ) mp <- sampleStanModel( - jm, - data = stan_data, - iter_sampling = 500, - iter_warmup = 1000, - chains = 1, - parallel_chains = 1, - exe_file = file.path("local", "full") + jm, + data = stan_data, + iter_sampling = 500, + iter_warmup = 1000, + chains = 1, + parallel_chains = 1, + exe_file = file.path("local", "full") ) diff --git a/misc/tests/gsf-with-link.R b/misc/tests/gsf-with-link.R index a7414b6d..bf8d16d5 100644 --- a/misc/tests/gsf-with-link.R +++ b/misc/tests/gsf-with-link.R @@ -18,36 +18,36 @@ devtools::load_all(export_all = FALSE) ## Generate Test data with known parameters jlist <- simulate_joint_data( - n_arm = c(80, 80), - times = seq(0, 4, by = (1/365)/2), - lambda_cen = 1 / 9000, - beta_cat = c( - "A" = 0, - "B" = -0.1, - "C" = 0.5 - ), - beta_cont = 0.3, - lm_fun = sim_lm_gsf( - sigma = 0.003, - mu_s = c(0.2, 0.25), - mu_g = c(0.15, 0.2), - mu_phi = c(0.4, 0.6), - mu_b = 60, - eta_b_sigma = 0.5, - eta_s_sigma = 0.5, - eta_g_sigma = 0.5, - eta_phi_sigma = 0.5, - omega_b = 0.1, - omega_s = 0.1, - omega_g = 0.1, - omega_phi = 0.2, - link_dsld = 0.02, - link_ttg = -0.02 - ), - os_fun = sim_os_weibull( - lambda = 365 * (1/400), - gamma = 1 - ) + n_arm = c(80, 80), + times = seq(0, 4, by = (1 / 365) / 2), + lambda_cen = 1 / 9000, + beta_cat = c( + "A" = 0, + "B" = -0.1, + "C" = 0.5 + ), + beta_cont = 0.3, + lm_fun = sim_lm_gsf( + sigma = 0.003, + mu_s = c(0.2, 0.25), + mu_g = c(0.15, 0.2), + mu_phi = c(0.4, 0.6), + mu_b = 60, + eta_b_sigma = 0.5, + eta_s_sigma = 0.5, + eta_g_sigma = 0.5, + eta_phi_sigma = 0.5, + omega_b = 0.1, + omega_s = 0.1, + omega_g = 0.1, + omega_phi = 0.2, + link_dsld = 0.02, + link_ttg = -0.02 + ), + os_fun = sim_os_weibull( + lambda = 365 * (1 / 400), + gamma = 1 + ) ) @@ -57,8 +57,8 @@ dat_os <- jlist$os select_times <- c(1, 100, 150, 200, 300, 400, 500, 600, 800, 900) * (1 / 365) dat_lm <- jlist$lm |> - dplyr::filter(time %in% select_times) |> - dplyr::arrange(pt, time) + dplyr::filter(time %in% select_times) |> + dplyr::arrange(pt, time) @@ -66,33 +66,29 @@ dat_lm <- jlist$lm |> pnam <- unique(dat_os$pt) |> sample(size = 10) ggplot(data = dat_lm |> dplyr::filter(pt %in% pnam)) + - geom_point(aes(x = time, y = sld, col =pt, group =pt)) + - geom_line(aes(x = time, y = sld, col =pt, group =pt)) + - theme_bw() + geom_point(aes(x = time, y = sld, col = pt, group = pt)) + + geom_line(aes(x = time, y = sld, col = pt, group = pt)) + + theme_bw() jm <- JointModel( - longitudinal = LongitudinalGSF( - - mu_bsld = Parameter(prior_lognormal(log(70), 5), init = 70), - mu_ks = Parameter(prior_lognormal(log(0.2), 1), init = 0.3), - mu_kg = Parameter(prior_lognormal( log(0.2), 1), init = 0.2), - mu_phi = Parameter(prior_beta(2, 2), init = 0.2), - - omega_bsld = Parameter(prior_lognormal(log(0.135), 1), init = 0.01), - omega_ks = Parameter(prior_lognormal(log(0.15), 1), init = 0.01), - omega_kg = Parameter(prior_lognormal(log(0.225), 1), init = 0.01), - omega_phi = Parameter(prior_lognormal(log(0.75), 1), init = 0.01), - - sigma = Parameter(prior_lognormal(log(0.01), 1), init = 0.01), - - tilde_bsld = Parameter(prior_normal(0, 5), init = 0.1), - tilde_ks = Parameter(prior_normal(0, 2), init = 0.1), - tilde_kg = Parameter(prior_normal(0, 1), init = 0.1), - tilde_phi = Parameter(prior_normal(0, 5), init = 0.1) - ), - survival = SurvivalExponential(), - link = LinkGSF() + longitudinal = LongitudinalGSF( + mu_bsld = Parameter(prior_lognormal(log(70), 5), init = 70), + mu_ks = Parameter(prior_lognormal(log(0.2), 1), init = 0.3), + mu_kg = Parameter(prior_lognormal(log(0.2), 1), init = 0.2), + mu_phi = Parameter(prior_beta(2, 2), init = 0.2), + omega_bsld = Parameter(prior_lognormal(log(0.135), 1), init = 0.01), + omega_ks = Parameter(prior_lognormal(log(0.15), 1), init = 0.01), + omega_kg = Parameter(prior_lognormal(log(0.225), 1), init = 0.01), + omega_phi = Parameter(prior_lognormal(log(0.75), 1), init = 0.01), + sigma = Parameter(prior_lognormal(log(0.01), 1), init = 0.01), + tilde_bsld = Parameter(prior_normal(0, 5), init = 0.1), + tilde_ks = Parameter(prior_normal(0, 2), init = 0.1), + tilde_kg = Parameter(prior_normal(0, 1), init = 0.1), + tilde_phi = Parameter(prior_normal(0, 5), init = 0.1) + ), + survival = SurvivalExponential(), + link = LinkGSF() ) @@ -105,31 +101,31 @@ write_stan(jm, "local/debug.stan") ## Prepare data for sampling jdat <- DataJoint( - survival = DataSurvival( - data = dat_os, - formula = Surv(time, event) ~ cov_cat + cov_cont, - subject = "pt", - arm = "arm", - study = "study" - ), - longitudinal = DataLongitudinal( - data = dat_lm, - formula = sld ~ time, - subject = "pt", - threshold = 5 - ) + survival = DataSurvival( + data = dat_os, + formula = Surv(time, event) ~ cov_cat + cov_cont, + subject = "pt", + arm = "arm", + study = "study" + ), + longitudinal = DataLongitudinal( + data = dat_lm, + formula = sld ~ time, + subject = "pt", + threshold = 5 + ) ) ## Sample from JointModel mp <- sampleStanModel( - jm, - data = jdat, - iter_sampling = 500, - iter_warmup = 1000, - chains = 1, - parallel_chains = 1, - exe_file = file.path("local", "full") + jm, + data = jdat, + iter_sampling = 500, + iter_warmup = 1000, + chains = 1, + parallel_chains = 1, + exe_file = file.path("local", "full") ) @@ -137,14 +133,14 @@ mp <- sampleStanModel( vars <- c( - "lm_gsf_mu_bsld[1]", - "lm_gsf_mu_phi[1]", "lm_gsf_mu_phi[2]", - "lm_gsf_mu_kg[1]", "lm_gsf_mu_kg[2]", - "lm_gsf_mu_ks[1]", "lm_gsf_mu_ks[2]", - "lm_gsf_sigma", - "lm_gsf_omega_bsld", "lm_gsf_omega_kg", - "lm_gsf_omega_phi", "lm_gsf_omega_ks", - "sm_exp_lambda", "lm_gsf_beta", "lm_gsf_gamma" + "lm_gsf_mu_bsld[1]", + "lm_gsf_mu_phi[1]", "lm_gsf_mu_phi[2]", + "lm_gsf_mu_kg[1]", "lm_gsf_mu_kg[2]", + "lm_gsf_mu_ks[1]", "lm_gsf_mu_ks[2]", + "lm_gsf_sigma", + "lm_gsf_omega_bsld", "lm_gsf_omega_kg", + "lm_gsf_omega_phi", "lm_gsf_omega_ks", + "sm_exp_lambda", "lm_gsf_beta", "lm_gsf_gamma" ) diff --git a/misc/tests/os-exponential.R b/misc/tests/os-exponential.R index 1699abc4..63beb50d 100644 --- a/misc/tests/os-exponential.R +++ b/misc/tests/os-exponential.R @@ -1,67 +1,66 @@ - devtools::document() devtools::load_all() library(bayesplot) library(survival) jlist <- simulate_joint_data( - n_arm = c(1000, 1000), - times = seq(1, 1000, by = 0.5), - lambda_cen = 1 / 9000, - beta_cat = c( - "A" = 0, - "B" = -0.3, - "C" = 0.5 - ), - beta_cont = 0.2, - lm_fun = sim_lm_random_slope(phi = 0), - os_fun = sim_os_exponential(lambda = 1/100) + n_arm = c(1000, 1000), + times = seq(1, 1000, by = 0.5), + lambda_cen = 1 / 9000, + beta_cat = c( + "A" = 0, + "B" = -0.3, + "C" = 0.5 + ), + beta_cont = 0.2, + lm_fun = sim_lm_random_slope(phi = 0), + os_fun = sim_os_exponential(lambda = 1 / 100) ) dat_os <- jlist$os dat_lm <- jlist$lm |> - dplyr::filter(time %in% c(1, 100, 150, 200, 300, 400, 500, 600, 800, 900)) |> - dplyr::arrange(time, pt) + dplyr::filter(time %in% c(1, 100, 150, 200, 300, 400, 500, 600, 800, 900)) |> + dplyr::arrange(time, pt) jm <- JointModel( - survival = SurvivalExponential() + survival = SurvivalExponential() ) write_stan(jm, "local/debug.stan") jdat <- DataJoint( - survival = DataSurvival( - data = dat_os, - formula = Surv(time, event) ~ cov_cat + cov_cont, - subject = "pt", - arm = "arm", - study = "study" - ), - longitudinal = DataLongitudinal( - data = dat_lm, - formula = sld ~ time, - subject = "pt", - threshold = 5 - ) + survival = DataSurvival( + data = dat_os, + formula = Surv(time, event) ~ cov_cat + cov_cont, + subject = "pt", + arm = "arm", + study = "study" + ), + longitudinal = DataLongitudinal( + data = dat_lm, + formula = sld ~ time, + subject = "pt", + threshold = 5 + ) ) mp <- sampleStanModel( - jm, - data = jdat, - iter_sampling = 500, - iter_warmup = 500, - chains = 1, - parallel_chains = 1, - exe_file = file.path("local", "full") + jm, + data = jdat, + iter_sampling = 500, + iter_warmup = 500, + chains = 1, + parallel_chains = 1, + exe_file = file.path("local", "full") ) vars <- c( - "sm_exp_lambda", # 0.01 - "beta_os_cov[1]", # -0.3 - "beta_os_cov[2]", # 0.5 - "beta_os_cov[3]" # 0.2 + "sm_exp_lambda", # 0.01 + "beta_os_cov[1]", # -0.3 + "beta_os_cov[2]", # 0.5 + "beta_os_cov[3]" # 0.2 ) mp$summary(vars) @@ -69,10 +68,3 @@ mp$summary(vars) mcmc_trace(mp$draws("sm_exp_lambda")) - - - - - - - diff --git a/misc/tests/os-loglogistic.R b/misc/tests/os-loglogistic.R index 5f3b1bc3..a6cdaabe 100644 --- a/misc/tests/os-loglogistic.R +++ b/misc/tests/os-loglogistic.R @@ -1,69 +1,68 @@ - devtools::document() devtools::load_all() library(bayesplot) library(survival) jlist <- simulate_joint_data( - n_arm = c(1000, 1000), - times = 1:2000, - lambda_cen = 1 / 9000, - beta_cat = c( - "A" = 0, - "B" = -0.1, - "C" = 0.5 - ), - beta_cont = 0.3, - lm_fun = sim_lm_random_slope(phi = 0), - os_fun = sim_os_loglogistic( - lambda = 1/50, - p = 1.1 - ) + n_arm = c(1000, 1000), + times = 1:2000, + lambda_cen = 1 / 9000, + beta_cat = c( + "A" = 0, + "B" = -0.1, + "C" = 0.5 + ), + beta_cont = 0.3, + lm_fun = sim_lm_random_slope(phi = 0), + os_fun = sim_os_loglogistic( + lambda = 1 / 50, + p = 1.1 + ) ) dat_os <- jlist$os dat_lm <- jlist$lm |> - dplyr::filter(time %in% c(1, 100, 150, 200, 300, 400, 500, 600, 800, 900)) |> - dplyr::arrange(time, pt) + dplyr::filter(time %in% c(1, 100, 150, 200, 300, 400, 500, 600, 800, 900)) |> + dplyr::arrange(time, pt) jm <- JointModel( - survival = SurvivalLogLogistic() + survival = SurvivalLogLogistic() ) write_stan(jm, "local/debug.stan") jdat <- DataJoint( - survival = DataSurvival( - data = dat_os, - formula = Surv(time, event) ~ cov_cat + cov_cont, - subject = "pt", - arm = "arm", - study = "study" - ), - longitudinal = DataLongitudinal( - data = dat_lm, - formula = sld ~ time, - subject = "pt", - threshold = 5 - ) + survival = DataSurvival( + data = dat_os, + formula = Surv(time, event) ~ cov_cat + cov_cont, + subject = "pt", + arm = "arm", + study = "study" + ), + longitudinal = DataLongitudinal( + data = dat_lm, + formula = sld ~ time, + subject = "pt", + threshold = 5 + ) ) mp <- sampleStanModel( - jm, - data = jdat, - iter_sampling = 1000, - iter_warmup = 500, - chains = 1, - exe_file = file.path("local", "full") + jm, + data = jdat, + iter_sampling = 1000, + iter_warmup = 500, + chains = 1, + exe_file = file.path("local", "full") ) vars <- c( - "sm_logl_lambda", # 0.02 - "sm_logl_p", # 1.1 - "beta_os_cov" # -0.1, 0.5, 0.3 + "sm_logl_lambda", # 0.02 + "sm_logl_p", # 1.1 + "beta_os_cov" # -0.1, 0.5, 0.3 ) mp$summary(vars) @@ -80,5 +79,3 @@ mcmc_trace(mp$draws("sm_logl_p")) # Surv(time, event) ~ cov_cat + cov_cont, # data = dat_os # ) - - diff --git a/misc/tests/os-weibull.R b/misc/tests/os-weibull.R index 16d564bf..dd0066ca 100644 --- a/misc/tests/os-weibull.R +++ b/misc/tests/os-weibull.R @@ -1,69 +1,68 @@ - devtools::document() devtools::load_all() library(bayesplot) library(survival) jlist <- simulate_joint_data( - n_arm = c(1000, 1000), - times = 1:2000, - lambda_cen = 1 / 9000, - beta_cat = c( - "A" = 0, - "B" = -0.1, - "C" = 0.5 - ), - beta_cont = 0.3, - lm_fun = sim_lm_random_slope(phi = 0), - os_fun = sim_os_weibull( - lambda = 0.00333, - gamma = 0.97 - ) + n_arm = c(1000, 1000), + times = 1:2000, + lambda_cen = 1 / 9000, + beta_cat = c( + "A" = 0, + "B" = -0.1, + "C" = 0.5 + ), + beta_cont = 0.3, + lm_fun = sim_lm_random_slope(phi = 0), + os_fun = sim_os_weibull( + lambda = 0.00333, + gamma = 0.97 + ) ) dat_os <- jlist$os dat_lm <- jlist$lm |> - dplyr::filter(time %in% c(1, 100, 150, 200, 300, 400, 500, 600, 800, 900)) |> - dplyr::arrange(time, pt) + dplyr::filter(time %in% c(1, 100, 150, 200, 300, 400, 500, 600, 800, 900)) |> + dplyr::arrange(time, pt) jm <- JointModel( - survival = SurvivalWeibullPH() + survival = SurvivalWeibullPH() ) write_stan(jm, "local/debug.stan") jdat <- DataJoint( - survival = DataSurvival( - data = dat_os, - formula = Surv(time, event) ~ cov_cat + cov_cont, - subject = "pt", - arm = "arm", - study = "study" - ), - longitudinal = DataLongitudinal( - data = dat_lm, - formula = sld ~ time, - subject = "pt", - threshold = 5 - ) + survival = DataSurvival( + data = dat_os, + formula = Surv(time, event) ~ cov_cat + cov_cont, + subject = "pt", + arm = "arm", + study = "study" + ), + longitudinal = DataLongitudinal( + data = dat_lm, + formula = sld ~ time, + subject = "pt", + threshold = 5 + ) ) mp <- sampleStanModel( - jm, - data = jdat, - iter_sampling = 1500, - iter_warmup = 1000, - chains = 1, - exe_file = file.path("local", "full") + jm, + data = jdat, + iter_sampling = 1500, + iter_warmup = 1000, + chains = 1, + exe_file = file.path("local", "full") ) vars <- c( - "sm_weibull_ph_lambda", # 0.00333 - "sm_weibull_ph_gamma", # 0.97 - "beta_os_cov" # -0.1, 0.5, 0.3 + "sm_weibull_ph_lambda", # 0.00333 + "sm_weibull_ph_gamma", # 0.97 + "beta_os_cov" # -0.1, 0.5, 0.3 ) mp$summary(vars) @@ -76,8 +75,6 @@ mcmc_trace(mp$draws("sm_weibull_ph_gamma")) library(survival) coxph( - Surv(time, event) ~ cov_cat + cov_cont, - data = dat_os + Surv(time, event) ~ cov_cat + cov_cont, + data = dat_os ) - - diff --git a/misc/tests/rs-no-link.R b/misc/tests/rs-no-link.R index f4154bab..3c5e6815 100644 --- a/misc/tests/rs-no-link.R +++ b/misc/tests/rs-no-link.R @@ -11,11 +11,11 @@ devtools::load_all(export_all = FALSE) #### Example 1 - Fully specified model jm <- JointModel( - longitudinal = LongitudinalRandomSlope( - intercept = prior_normal(30, 2), - slope_sigma = prior_lognormal(0.2, sigma = 3), - sigma = prior_normal(3, sigma = 0.5) - ) + longitudinal = LongitudinalRandomSlope( + intercept = prior_normal(30, 2), + slope_sigma = prior_lognormal(0.2, sigma = 3), + sigma = prior_normal(3, sigma = 0.5) + ) ) @@ -26,73 +26,73 @@ write_stan(jm, "local/debug.stan") ## Generate Test data with known parameters jlist <- simulate_joint_data( - n = c(400, 400), - times = 1:2000, - lambda_cen = 1 / 9000, - beta_cat = c( - "A" = 0, - "B" = -0.1, - "C" = 0.5 - ), - beta_cont = 0.3, - lm_fun = sim_lm_random_slope( - intercept = 30, - sigma = 3, - slope_mu = c(1,3), - slope_sigma = 0.2, - phi = 0, - .debug = TRUE - ), - os_fun = sim_os_exponential( - lambda = 0.00333 # 1/300 - ) + n = c(400, 400), + times = 1:2000, + lambda_cen = 1 / 9000, + beta_cat = c( + "A" = 0, + "B" = -0.1, + "C" = 0.5 + ), + beta_cont = 0.3, + lm_fun = sim_lm_random_slope( + intercept = 30, + sigma = 3, + slope_mu = c(1, 3), + slope_sigma = 0.2, + phi = 0, + .debug = TRUE + ), + os_fun = sim_os_exponential( + lambda = 0.00333 # 1/300 + ) ) ## Extract data to individual datasets dat_os <- jlist$os dat_lm <- jlist$lm |> - dplyr::filter(time %in% c(1, 50, 100, 150, 200, 250, 300)) |> - dplyr::arrange(time, pt) + dplyr::filter(time %in% c(1, 50, 100, 150, 200, 250, 300)) |> + dplyr::arrange(time, pt) ## Prepare data for sampling jdat <- DataJoint( - survival = DataSurvival( - data = dat_os, - formula = Surv(time, event) ~ cov_cat + cov_cont, - subject = "pt", - arm = "arm", - study = "study" - ), - longitudinal = DataLongitudinal( - data = dat_lm, - formula = sld ~ time, - subject = "pt", - threshold = 5 - ) + survival = DataSurvival( + data = dat_os, + formula = Surv(time, event) ~ cov_cat + cov_cont, + subject = "pt", + arm = "arm", + study = "study" + ), + longitudinal = DataLongitudinal( + data = dat_lm, + formula = sld ~ time, + subject = "pt", + threshold = 5 + ) ) ## Sample from JointModel mp <- sampleStanModel( - jm, - data = jdat, - iter_sampling = 1000, - iter_warmup = 1000, - chains = 1, - parallel_chains = 1, - exe_file = file.path("local", "full") + jm, + data = jdat, + iter_sampling = 1000, + iter_warmup = 1000, + chains = 1, + parallel_chains = 1, + exe_file = file.path("local", "full") ) ### Select which parameters we actually care about ### Not all will exist depending on which model was run vars <- c( - "lm_rs_intercept", # 30 - "lm_rs_slope_mu", # 1 , 3 - "lm_rs_slope_sigma", # 0.2 - "lm_rs_sigma" # 3 + "lm_rs_intercept", # 30 + "lm_rs_slope_mu", # 1 , 3 + "lm_rs_slope_sigma", # 0.2 + "lm_rs_sigma" # 3 ) mp$summary(vars) @@ -105,16 +105,16 @@ mp$summary(vars) # slopes <- tibble( - est = mp$summary("lm_rs_rslope")$mean, - actual = dat_lm |> group_by(pt) |> slice(1) |> pull(slope_ind) + est = mp$summary("lm_rs_rslope")$mean, + actual = dat_lm |> group_by(pt) |> slice(1) |> pull(slope_ind) ) mean(slopes$actual - slopes$est) sd(slopes$actual - slopes$est) ggplot(slopes, aes(x = est, y = actual)) + - geom_point() + - geom_abline() + geom_point() + + geom_abline() @@ -130,17 +130,7 @@ coxph(data = dat_os, Surv(time, event) ~ cov_cat + cov_cont) mod <- survreg(Surv(time, event) ~ cov_cat + cov_cont, data = dat_os, dist = "weibull") c( - "lambda" = exp(-coef(mod)[["(Intercept)"]]/mod$scale), - "gamma" = 1/mod$scale, - -coef(mod)[-1]/mod$scale + "lambda" = exp(-coef(mod)[["(Intercept)"]] / mod$scale), + "gamma" = 1 / mod$scale, + -coef(mod)[-1] / mod$scale ) - - - - - - - - - - diff --git a/misc/tests/rs-with-link.R b/misc/tests/rs-with-link.R index a76f7a54..f540f36f 100644 --- a/misc/tests/rs-with-link.R +++ b/misc/tests/rs-with-link.R @@ -11,9 +11,9 @@ devtools::load_all(export_all = FALSE) #### Example 1 - Fully specified model - using the defaults for everything jm <- JointModel( - longitudinal = LongitudinalRandomSlope(), - survival = SurvivalExponential(), - link = LinkRandomSlope() + longitudinal = LongitudinalRandomSlope(), + survival = SurvivalExponential(), + link = LinkRandomSlope() ) @@ -24,76 +24,76 @@ write_stan(jm, "local/debug.stan") ## Generate Test data with known parameters jlist <- simulate_joint_data( - n = c(400, 400), - times = 1:2000, - lambda_cen = 1 / 9000, - beta_cat = c( - "A" = 0, - "B" = -0.1, - "C" = 0.5 - ), - beta_cont = 0.3, - lm_fun = sim_lm_random_slope( - intercept = 30, - sigma = 3, - slope_mu = c(1,3), - slope_sigma = 0.2, - phi = 0.1, - .debug = TRUE - ), - os_fun = sim_os_exponential( - lambda = 0.00333 # 1/300 - ) + n = c(400, 400), + times = 1:2000, + lambda_cen = 1 / 9000, + beta_cat = c( + "A" = 0, + "B" = -0.1, + "C" = 0.5 + ), + beta_cont = 0.3, + lm_fun = sim_lm_random_slope( + intercept = 30, + sigma = 3, + slope_mu = c(1, 3), + slope_sigma = 0.2, + phi = 0.1, + .debug = TRUE + ), + os_fun = sim_os_exponential( + lambda = 0.00333 # 1/300 + ) ) ## Extract data to individual datasets dat_os <- jlist$os dat_lm <- jlist$lm |> - dplyr::filter(time %in% c(1, 50, 100, 150, 200, 250, 300)) |> - dplyr::arrange(time, pt) + dplyr::filter(time %in% c(1, 50, 100, 150, 200, 250, 300)) |> + dplyr::arrange(time, pt) ## Prepare data for sampling jdat <- DataJoint( - survival = DataSurvival( - data = dat_os, - formula = Surv(time, event) ~ cov_cat + cov_cont, - subject = "pt", - arm = "arm", - study = "study" - ), - longitudinal = DataLongitudinal( - data = dat_lm, - formula = sld ~ time, - subject = "pt", - threshold = 5 - ) + survival = DataSurvival( + data = dat_os, + formula = Surv(time, event) ~ cov_cat + cov_cont, + subject = "pt", + arm = "arm", + study = "study" + ), + longitudinal = DataLongitudinal( + data = dat_lm, + formula = sld ~ time, + subject = "pt", + threshold = 5 + ) ) ## Sample from JointModel mp <- sampleStanModel( - jm, - data = jdat, - iter_sampling = 1000, - iter_warmup = 1000, - chains = 1, - parallel_chains = 1, - exe_file = file.path("local", "full") + jm, + data = jdat, + iter_sampling = 1000, + iter_warmup = 1000, + chains = 1, + parallel_chains = 1, + exe_file = file.path("local", "full") ) ### Select which parameters we actually care about ### Not all will exist depending on which model was run vars <- c( - "sm_exp_lambda", # 0.00333 - "beta_os_cov", # -0.1, 0.5, 0.3 - "lm_rs_intercept", # 30 - "lm_rs_slope_mu", # 1 , 3 - "lm_rs_slope_sigma", # 0.2 - "lm_rs_sigma", # 3 - "link_lm_phi" # 0.1 + "sm_exp_lambda", # 0.00333 + "beta_os_cov", # -0.1, 0.5, 0.3 + "lm_rs_intercept", # 30 + "lm_rs_slope_mu", # 1 , 3 + "lm_rs_slope_sigma", # 0.2 + "lm_rs_sigma", # 3 + "link_lm_phi" # 0.1 ) mp$summary(vars) @@ -106,16 +106,16 @@ mp$summary(vars) # slopes <- tibble( - est = mp$summary("lm_rs_rslope")$mean, - actual = dat_lm |> group_by(pt) |> slice(1) |> pull(slope_ind) + est = mp$summary("lm_rs_rslope")$mean, + actual = dat_lm |> group_by(pt) |> slice(1) |> pull(slope_ind) ) mean(slopes$actual - slopes$est) sd(slopes$actual - slopes$est) ggplot(slopes, aes(x = est, y = actual)) + - geom_point() + - geom_abline() + geom_point() + + geom_abline() @@ -131,17 +131,7 @@ coxph(data = dat_os, Surv(time, event) ~ cov_cat + cov_cont) mod <- survreg(Surv(time, event) ~ cov_cat + cov_cont, data = dat_os, dist = "weibull") c( - "lambda" = exp(-coef(mod)[["(Intercept)"]]/mod$scale), - "gamma" = 1/mod$scale, - -coef(mod)[-1]/mod$scale + "lambda" = exp(-coef(mod)[["(Intercept)"]] / mod$scale), + "gamma" = 1 / mod$scale, + -coef(mod)[-1] / mod$scale ) - - - - - - - - - - diff --git a/misc/tests/stan_functions.R b/misc/tests/stan_functions.R index aef4de29..fd2bee02 100644 --- a/misc/tests/stan_functions.R +++ b/misc/tests/stan_functions.R @@ -20,33 +20,33 @@ model { " run_stan_function <- function(stan_insert, stan_data, fun_files) { - stan_mods <- lapply(fun_files, \(i) StanModule(i)) - stan_mod <- Reduce(merge, stan_mods) - stan_plus_user <- merge( - stan_mod, - StanModule(stan_insert) - ) - stan_data$y_mean <- 0 - stan_final <- merge( - stan_plus_user, - StanModule(NULL_MODEL_INSERT) - ) - mod <- cmdstan_model(write_stan_file(as.character(stan_final))) - devnull <- capture.output({ - suppressMessages({ - suppressWarnings({ - fit <- mod$sample( - data = stan_data, - refresh = 0, - chains = 1, - iter_sampling = 1, - iter_warmup = 1, - show_messages = FALSE - ) - }) - }) + stan_mods <- lapply(fun_files, \(i) StanModule(i)) + stan_mod <- Reduce(merge, stan_mods) + stan_plus_user <- merge( + stan_mod, + StanModule(stan_insert) + ) + stan_data$y_mean <- 0 + stan_final <- merge( + stan_plus_user, + StanModule(NULL_MODEL_INSERT) + ) + mod <- cmdstan_model(write_stan_file(as.character(stan_final))) + devnull <- capture.output({ + suppressMessages({ + suppressWarnings({ + fit <- mod$sample( + data = stan_data, + refresh = 0, + chains = 1, + iter_sampling = 1, + iter_warmup = 1, + show_messages = FALSE + ) + }) }) - return(fit) + }) + return(fit) } @@ -59,12 +59,12 @@ data { vector[sld_n] sld_s; vector[sld_n] sld_g; vector[sld_n] sld_phi; - + int nrm_logden_n; vector[nrm_logden_n] nrm_logden_obs; vector[nrm_logden_n] nrm_logden_mu; vector[nrm_logden_n] nrm_logden_sigma; - + int nrm_logcum_n; real nrm_logcum_quant; vector[nrm_logcum_n] nrm_logcum_mu; @@ -74,13 +74,13 @@ data { generated quantities { vector[sld_n] sld_results = sld(sld_time, sld_bsld, sld_s, sld_g, sld_phi); - + vector[nrm_logden_n] nrm_logden_results = vect_normal_log_dens( nrm_logden_obs, nrm_logden_mu, nrm_logden_sigma ); - + vector[nrm_logcum_n] nrm_logcum_results = vect_normal_log_cum( nrm_logcum_quant, nrm_logcum_mu, @@ -90,19 +90,17 @@ generated quantities { " stan_data <- list( - sld_time = c(0, 0.5, 1, 2), - sld_bsld = c(50, 40, 40, 50), - sld_s = c(0.1, 0.2, 0.3, 0.4), - sld_g = c(0.2, 0.3, 0.2, 0.3), - sld_phi = c(0.5, 0.4, 0.5, 0.3), - - nrm_logden_obs = c(5, 4, 3, 9, 10), - nrm_logden_mu = c(4, 4, 3, 3, 10), - nrm_logden_sigma = c(2, 3, 4, 5, 6), - - nrm_logcum_quant = 5, - nrm_logcum_mu = c(4, 4, 3, 3, 10), - nrm_logcum_sigma = c(2, 3, 4, 5, 6) + sld_time = c(0, 0.5, 1, 2), + sld_bsld = c(50, 40, 40, 50), + sld_s = c(0.1, 0.2, 0.3, 0.4), + sld_g = c(0.2, 0.3, 0.2, 0.3), + sld_phi = c(0.5, 0.4, 0.5, 0.3), + nrm_logden_obs = c(5, 4, 3, 9, 10), + nrm_logden_mu = c(4, 4, 3, 3, 10), + nrm_logden_sigma = c(2, 3, 4, 5, 6), + nrm_logcum_quant = 5, + nrm_logcum_mu = c(4, 4, 3, 3, 10), + nrm_logcum_sigma = c(2, 3, 4, 5, 6) ) stan_data["sld_n"] <- length(stan_data$sld_time) @@ -110,9 +108,9 @@ stan_data["nrm_logden_n"] <- length(stan_data$nrm_logden_mu) stan_data["nrm_logcum_n"] <- length(stan_data$nrm_logcum_mu) fit <- run_stan_function( - stan_insert, - stan_data, - c("base/functions.stan", "lm-gsf/functions.stan") + stan_insert, + stan_data, + c("base/functions.stan", "lm-gsf/functions.stan") ) @@ -123,13 +121,13 @@ fit <- run_stan_function( observed <- as.numeric(fit$draws(variables = "sld_results")) |> round(4) expected <- sld( - stan_data$sld_time, - stan_data$sld_bsld, - stan_data$sld_s, - stan_data$sld_g, - stan_data$sld_phi + stan_data$sld_time, + stan_data$sld_bsld, + stan_data$sld_s, + stan_data$sld_g, + stan_data$sld_phi ) |> - round(4) + round(4) expect_equal(observed, expected) @@ -142,12 +140,12 @@ expect_equal(observed, expected) observed <- as.numeric(fit$draws(variables = "nrm_logden_results")) |> round(5) expected <- dnorm( - stan_data$nrm_logden_obs, - stan_data$nrm_logden_mu, - stan_data$nrm_logden_sigma, - log = TRUE + stan_data$nrm_logden_obs, + stan_data$nrm_logden_mu, + stan_data$nrm_logden_sigma, + log = TRUE ) |> - round(5) + round(5) expect_equal(observed, expected) @@ -159,12 +157,11 @@ expect_equal(observed, expected) observed <- as.numeric(fit$draws(variables = "nrm_logcum_results")) |> round(5) expected <- pnorm( - stan_data$nrm_logcum_quant, - stan_data$nrm_logcum_mu, - stan_data$nrm_logcum_sigma, - log = TRUE + stan_data$nrm_logcum_quant, + stan_data$nrm_logcum_mu, + stan_data$nrm_logcum_sigma, + log = TRUE ) |> - round(5) + round(5) expect_equal(observed, expected) - diff --git a/tests/testthat/test-DataJoint.R b/tests/testthat/test-DataJoint.R index 54b25009..8f5bc51c 100644 --- a/tests/testthat/test-DataJoint.R +++ b/tests/testthat/test-DataJoint.R @@ -2,100 +2,88 @@ library(survival) test_that("DataJoint errors if subjects don't allign after", { - df_surv <- data.frame( - vpt = c("A", "B", "C"), - vtime = c(100, 200, 150), - vevent = c(0, 1, 1), - vct = c(5, 2, 4), - varm = c("A1", "A1", "A2"), - vstudy = c("S1", "S1", "S2") + df_surv <- data.frame( + vpt = c("A", "B", "C"), + vtime = c(100, 200, 150), + vevent = c(0, 1, 1), + vct = c(5, 2, 4), + varm = c("A1", "A1", "A2"), + vstudy = c("S1", "S1", "S2") + ) + + df_long <- data.frame( + vpt = c("A", "A", "B", "B", "C", "B"), + vtime = c(10, 10, 20, 30, 40, 50), + vout = c(1, 2, 3, 4, 5, 6) + ) + + + do_surv <- DataSurvival( + data = df_surv, + formula = Surv(vtime, vevent) ~ vct, + subject = "vpt", + arm = "varm", + study = "vstudy" + ) + + do_long <- DataLongitudinal( + data = df_long, + formula = vout ~ vtime, + subject = "vpt" + ) + + d_joint <- DataJoint( + survival = do_surv, + longitudinal = do_long + ) + + ## Smoke tests + expect_equal(as.list(d_joint)$Nind, 3) + expect_equal(as.list(d_joint)$Nta_total, 6) + + + + df_long2 <- df_long + df_long2$vpt[df_long2$vpt == "C"] <- NA_character_ + + df_surv2 <- df_surv + df_surv2$vpt[df_surv2$vpt == "C"] <- NA_character_ + + suppressMessages({ + do_long2 <- DataLongitudinal( + data = df_long2, + formula = vout ~ vtime, + subject = "vpt" ) - - df_long <- data.frame( - vpt = c("A", "A", "B", "B", "C", "B"), - vtime = c(10, 10, 20, 30, 40, 50), - vout = c(1,2,3,4,5, 6) - ) - - - do_surv <- DataSurvival( - data = df_surv, - formula = Surv(vtime, vevent) ~ vct, - subject = "vpt", - arm = "varm", - study = "vstudy" - ) - - do_long <- DataLongitudinal( - data = df_long, - formula = vout ~ vtime, - subject = "vpt" + }) + + suppressMessages({ + do_surv2 <- DataSurvival( + data = df_surv2, + formula = Surv(vtime, vevent) ~ vct, + subject = "vpt", + arm = "varm", + study = "vstudy" ) + }) - d_joint <- DataJoint( + expect_error( + { + d_joint <- DataJoint( survival = do_surv, + longitudinal = do_long2 + ) + }, + regexp = "subjects in the survival" + ) + + expect_error( + { + d_joint <- DataJoint( + survival = do_surv2, longitudinal = do_long - ) - - ## Smoke tests - expect_equal(as.list(d_joint)$Nind, 3) - expect_equal(as.list(d_joint)$Nta_total, 6) - - - - df_long2 <- df_long - df_long2$vpt[df_long2$vpt == "C"] <- NA_character_ - - df_surv2 <- df_surv - df_surv2$vpt[df_surv2$vpt == "C"] <- NA_character_ - - suppressMessages({ - do_long2 <- DataLongitudinal( - data = df_long2, - formula = vout ~ vtime, - subject = "vpt" - ) - }) - - suppressMessages({ - do_surv2 <- DataSurvival( - data = df_surv2, - formula = Surv(vtime, vevent) ~ vct, - subject = "vpt", - arm = "varm", - study = "vstudy" - ) - }) - - expect_error( - { - d_joint <- DataJoint( - survival = do_surv, - longitudinal = do_long2 - ) - }, - regexp = "subjects in the survival" - ) - - expect_error( - { - d_joint <- DataJoint( - survival = do_surv2, - longitudinal = do_long - ) - }, - regexp = "subjects in the longitudinal" - ) + ) + }, + regexp = "subjects in the longitudinal" + ) }) - - - - - - - - - - - - diff --git a/tests/testthat/test-DataLongitudinal.R b/tests/testthat/test-DataLongitudinal.R index 3bd86092..855d0015 100644 --- a/tests/testthat/test-DataLongitudinal.R +++ b/tests/testthat/test-DataLongitudinal.R @@ -1,31 +1,27 @@ - - - test_that("DataLongitudinal being rendered to list is as expected for simple inputs", { - x <- data.frame( - vpt = c("b", "a", "a", "b", "c", "a"), - vtime = c(10, 20, 15, 5, 25, 35), - voutcome = c(2, 1, 4, 5, 6, 2) - ) - - - - dl <- DataLongitudinal( - data = x, - formula = voutcome ~ vtime, - subject = "vpt", - threshold = 3 - ) - - li <- as.list(dl) - - expect_equal(li$Yobs, x$voutcome) - expect_equal(li$ind_index, c(2, 1, 1, 2, 3, 1)) - expect_equal(li$obs_y_index, c(3, 4, 5)) - expect_equal(li$cens_y_index, c(1, 2, 6)) - expect_equal(li$Nta_cens_y, 3) - expect_equal(li$Nta_obs_y, 3) - expect_equal(li$Nta_total, 6) - expect_equal(li$Tobs, x$vtime) + x <- data.frame( + vpt = c("b", "a", "a", "b", "c", "a"), + vtime = c(10, 20, 15, 5, 25, 35), + voutcome = c(2, 1, 4, 5, 6, 2) + ) + + + + dl <- DataLongitudinal( + data = x, + formula = voutcome ~ vtime, + subject = "vpt", + threshold = 3 + ) + + li <- as.list(dl) + + expect_equal(li$Yobs, x$voutcome) + expect_equal(li$ind_index, c(2, 1, 1, 2, 3, 1)) + expect_equal(li$obs_y_index, c(3, 4, 5)) + expect_equal(li$cens_y_index, c(1, 2, 6)) + expect_equal(li$Nta_cens_y, 3) + expect_equal(li$Nta_obs_y, 3) + expect_equal(li$Nta_total, 6) + expect_equal(li$Tobs, x$vtime) }) - diff --git a/tests/testthat/test-DataSurvival.R b/tests/testthat/test-DataSurvival.R index a19160ba..b8c19ca2 100644 --- a/tests/testthat/test-DataSurvival.R +++ b/tests/testthat/test-DataSurvival.R @@ -1,66 +1,60 @@ - library(survival) x <- data.frame( - vpt = c("b" , "a", "c"), - varm = c("a", "a", "a"), - vstudy = c("s1", "s1", "s1"), - vtime = c(10 , 20 , 30), - vevent = c(1,1,1) + vpt = c("b", "a", "c"), + varm = c("a", "a", "a"), + vstudy = c("s1", "s1", "s1"), + vtime = c(10, 20, 30), + vevent = c(1, 1, 1) ) test_that("Error Handling", { + x2 <- x + x2$vpt <- "b" - x2 <- x - x2$vpt <- "b" - - expect_error( - { - DataSurvival( - data = x2, - formula = Surv(vtime, vevent) ~ 1, - subject = "vpt", - arm = "varm", - study = "vstudy" - ) - }, - "Only 1 survival observation" - ) -}) - - - -test_that("as.list returns expected values", { - - x <- data.frame( - vpt = c("b" , "a", "c", "d", "e"), - varm = c("b", "a", "a", "b", "b"), - vstudy = c("s1", "s1", "s1", "s1", "s1"), - vtime = c(10 , 20 , 30, 25, 15), - vevent = c(1,1,0, 1, 0 ) - ) - - df <- DataSurvival( - data = x, + expect_error( + { + DataSurvival( + data = x2, formula = Surv(vtime, vevent) ~ 1, subject = "vpt", arm = "varm", study = "vstudy" - ) + ) + }, + "Only 1 survival observation" + ) +}) - res <- as.list(df) - expect_equal(res$Nind_dead, 3) - expect_equal(res$dead_ind_index, c(1,2, 4)) - expect_equal(res$n_arms, 2) - expect_equal(res$Nind, 5) - expect_equal(res$Times, x$vtime) - expect_equal(res$study_index, rep(1, 5)) - expect_equal(res$arm_index, c(2,1,1,2,2)) - expect_equal(res$n_index_per_arm, c(2, 3)) - expect_equal(res$index_per_arm, c(2,3, 1, 4, 5)) +test_that("as.list returns expected values", { + x <- data.frame( + vpt = c("b", "a", "c", "d", "e"), + varm = c("b", "a", "a", "b", "b"), + vstudy = c("s1", "s1", "s1", "s1", "s1"), + vtime = c(10, 20, 30, 25, 15), + vevent = c(1, 1, 0, 1, 0) + ) + + df <- DataSurvival( + data = x, + formula = Surv(vtime, vevent) ~ 1, + subject = "vpt", + arm = "varm", + study = "vstudy" + ) + + res <- as.list(df) + + expect_equal(res$Nind_dead, 3) + expect_equal(res$dead_ind_index, c(1, 2, 4)) + expect_equal(res$n_arms, 2) + expect_equal(res$Nind, 5) + expect_equal(res$Times, x$vtime) + expect_equal(res$study_index, rep(1, 5)) + expect_equal(res$arm_index, c(2, 1, 1, 2, 2)) + expect_equal(res$n_index_per_arm, c(2, 3)) + expect_equal(res$index_per_arm, c(2, 3, 1, 4, 5)) }) - - diff --git a/tests/testthat/test-JointModel.R b/tests/testthat/test-JointModel.R index 3b41fdba..47d17e12 100644 --- a/tests/testthat/test-JointModel.R +++ b/tests/testthat/test-JointModel.R @@ -1,21 +1,11 @@ - - - - test_that("JointModel smoke tests", { - jm <- JointModel( - longitudinal = LongitudinalRandomSlope(), - survival = SurvivalWeibullPH(), - link = LinkRandomSlope() - ) - - jm_char <- as.character(jm) - expect_equal(length(jm_char), 1) - expect_true(nchar(jm_char) > 3000) + jm <- JointModel( + longitudinal = LongitudinalRandomSlope(), + survival = SurvivalWeibullPH(), + link = LinkRandomSlope() + ) + + jm_char <- as.character(jm) + expect_equal(length(jm_char), 1) + expect_true(nchar(jm_char) > 3000) }) - - - - - - diff --git a/tests/testthat/test-LinkNone.R b/tests/testthat/test-LinkNone.R index f4d540e9..8a4e58ef 100644 --- a/tests/testthat/test-LinkNone.R +++ b/tests/testthat/test-LinkNone.R @@ -1,20 +1,13 @@ - - - test_that("LinkNone object is generated without unexpected input", { - expect_error( - LinkNone(stan = StanModule(), parameters = ParameterList(mu = 0.2)), - "unused argument", - ignore.case = TRUE - ) - - link <- LinkNone() - expect_true(is(link, "Link")) - expect_true(is(link, "LinkNone")) - expect_false(is(link, "LongitudinalModel")) - expect_false(is(link, "link_gsf_abstract")) + expect_error( + LinkNone(stan = StanModule(), parameters = ParameterList(mu = 0.2)), + "unused argument", + ignore.case = TRUE + ) + + link <- LinkNone() + expect_true(is(link, "Link")) + expect_true(is(link, "LinkNone")) + expect_false(is(link, "LongitudinalModel")) + expect_false(is(link, "link_gsf_abstract")) }) - - - - diff --git a/tests/testthat/test-LinkRandomSlope.R b/tests/testthat/test-LinkRandomSlope.R index a1834d06..3a16883e 100644 --- a/tests/testthat/test-LinkRandomSlope.R +++ b/tests/testthat/test-LinkRandomSlope.R @@ -1,34 +1,22 @@ - - - test_that("LinkRandomSlope smoke tests", { - - linkobject <- LinkRandomSlope( - link_lm_phi = prior_normal(0, 5, init = 0.01) - ) - - - expect_equal( - initialValues(linkobject), - list("link_lm_phi" = 0.01) - ) - - - expect_true( - is(as.StanModule(linkobject), "StanModule") - ) - - - expect_error( - LongitudinalRandomSlope(lm_rs_pp = prior_normal(0, 1)), - "unused argument \\(lm_rs_pp" - ) -}) - - + linkobject <- LinkRandomSlope( + link_lm_phi = prior_normal(0, 5, init = 0.01) + ) + expect_equal( + initialValues(linkobject), + list("link_lm_phi" = 0.01) + ) + expect_true( + is(as.StanModule(linkobject), "StanModule") + ) + expect_error( + LongitudinalRandomSlope(lm_rs_pp = prior_normal(0, 1)), + "unused argument \\(lm_rs_pp" + ) +}) diff --git a/tests/testthat/test-ParameterList.R b/tests/testthat/test-ParameterList.R index 4d207000..151a5272 100644 --- a/tests/testthat/test-ParameterList.R +++ b/tests/testthat/test-ParameterList.R @@ -1,32 +1,28 @@ +test_that("ParameterList smoke tests", { + pl <- ParameterList( + Parameter(name = "inter", prior = prior_gamma(1, 2)), + Parameter(name = "myp", prior = prior_normal(1, 4, init = 8)) + ) + # Can extract parameter names + expect_equal(names(pl), c("inter", "myp")) + # Can extract initial values + actual <- initialValues(pl) + expected <- list( + "inter" = 1 / 2, + "myp" = 8 + ) + expect_equal(actual, expected) -test_that("ParameterList smoke tests", { - pl <- ParameterList( - Parameter(name = "inter", prior = prior_gamma(1, 2)), - Parameter(name = "myp", prior = prior_normal(1, 4, init = 8)) - ) - - # Can extract parameter names - expect_equal(names(pl), c("inter", "myp")) - - - # Can extract initial values - actual <- initialValues(pl) - expected <- list( - "inter" = 1 / 2, - "myp" = 8 - ) - expect_equal(actual, expected) - - # Can render to character - actual <- c( - " inter ~ gamma(1, 2);", - " myp ~ normal(1, 4);" - ) |> - paste(collapse = "\n") + # Can render to character + actual <- c( + " inter ~ gamma(1, 2);", + " myp ~ normal(1, 4);" + ) |> + paste(collapse = "\n") - expect_equal(actual, as.character(pl)) + expect_equal(actual, as.character(pl)) }) diff --git a/tests/testthat/test-Parameters.R b/tests/testthat/test-Parameters.R index e8fdd32a..1411f440 100644 --- a/tests/testthat/test-Parameters.R +++ b/tests/testthat/test-Parameters.R @@ -1,22 +1,9 @@ - - - - - - - test_that("Parameters smoke tests", { - p <- Parameter(name = "intercept", prior = prior_beta(5, 4)) - expect_equal(initialValues(p), 5 / (5 + 4)) - expect_equal(names(p), "intercept") - - p <- Parameter(name = "myp", prior = prior_beta(5, 4, init = 9)) - expect_equal(initialValues(p), 9) - expect_equal(names(p), "myp") -}) - - - - - + p <- Parameter(name = "intercept", prior = prior_beta(5, 4)) + expect_equal(initialValues(p), 5 / (5 + 4)) + expect_equal(names(p), "intercept") + p <- Parameter(name = "myp", prior = prior_beta(5, 4, init = 9)) + expect_equal(initialValues(p), 9) + expect_equal(names(p), "myp") +}) diff --git a/tests/testthat/test-Prior.R b/tests/testthat/test-Prior.R index 98602e81..f60c5fbe 100644 --- a/tests/testthat/test-Prior.R +++ b/tests/testthat/test-Prior.R @@ -1,15 +1,9 @@ - - - test_that("Priors work as expected", { - x <- prior_normal(4, 10) - expect_equal(initialValues(x), 4) - expect_equal(as.character(x), "normal(4, 10);") - - x <- prior_normal(4, 10, 20) - expect_equal(initialValues(x), 20) - expect_equal(as.character(x), "normal(4, 10);") -}) - - + x <- prior_normal(4, 10) + expect_equal(initialValues(x), 4) + expect_equal(as.character(x), "normal(4, 10);") + x <- prior_normal(4, 10, 20) + expect_equal(initialValues(x), 20) + expect_equal(as.character(x), "normal(4, 10);") +}) diff --git a/tests/testthat/test-compile.R b/tests/testthat/test-compile.R index 2fa910b6..96b8130b 100644 --- a/tests/testthat/test-compile.R +++ b/tests/testthat/test-compile.R @@ -1,4 +1,3 @@ - model <- " data { int n; @@ -17,16 +16,9 @@ model { test_that("compileStanModel doesn't error if the directory doesn't exist", { + smod <- StanModule(model) + fpath <- file.path(tempdir(), "abcd", "efg", "model") + z <- compileStanModel(smod, fpath) - - - smod <- StanModule(model) - fpath <- file.path(tempdir(), "abcd", "efg", "model") - z <- compileStanModel(smod, fpath) - - expect_true(file.exists(fpath)) + expect_true(file.exists(fpath)) }) - - - - diff --git a/tests/testthat/test-utilities.R b/tests/testthat/test-utilities.R index 60aa24a8..2a2049f8 100644 --- a/tests/testthat/test-utilities.R +++ b/tests/testthat/test-utilities.R @@ -1,82 +1,72 @@ - - library(survival) test_that("get_missing_rownumbers works as expected", { - df <- data.frame( - time = c(1, 2, 3, NA), - event = c(0, NA, 1, 1), - x = c(NA, NA, NA, NA), - y = c(NA, 5, 3, 7), - z = c(1, 4, 3, NA) - ) - - actual <- get_missing_rownumbers(df, Surv(time, event) ~ y + z) - expected <- c(1, 2, 4) - expect_equal(actual, expected) + df <- data.frame( + time = c(1, 2, 3, NA), + event = c(0, NA, 1, 1), + x = c(NA, NA, NA, NA), + y = c(NA, 5, 3, 7), + z = c(1, 4, 3, NA) + ) + + actual <- get_missing_rownumbers(df, Surv(time, event) ~ y + z) + expected <- c(1, 2, 4) + expect_equal(actual, expected) }) test_that("remove_missing_rows works as expected", { - - df <- data.frame( - time = c(1 , NA, 3, 5, 5, 4, 5, 5), - event = c(0 , 0, NA, 1, 0, 1, 1, 0), - y = c(NA, 5, 3, NA, 3, 4, 3, 3), - z = c(NA, 4, 3, 2, NA, 2, 4, 3), - x = c(NA, NA, NA, NA, NA, NA, NA, NA), - w = c("a", "b", "c", "d", "f", "h", "i", NA_character_) - ) - - expected <- data.frame( - time = c( 4, 5), - event = c( 1, 1), - y = c( 4, 3), - z = c( 2, 4), - x = c(NA, NA), - w = c("h", "i") - ) - suppressMessages({ - actual <- remove_missing_rows(df, Surv(time, event) ~ y, extra_vars = c("z", "w")) - }) - - rownames(expected) <- NULL - rownames(actual) <- NULL - - expect_equal(expected, actual) - + df <- data.frame( + time = c(1, NA, 3, 5, 5, 4, 5, 5), + event = c(0, 0, NA, 1, 0, 1, 1, 0), + y = c(NA, 5, 3, NA, 3, 4, 3, 3), + z = c(NA, 4, 3, 2, NA, 2, 4, 3), + x = c(NA, NA, NA, NA, NA, NA, NA, NA), + w = c("a", "b", "c", "d", "f", "h", "i", NA_character_) + ) + + expected <- data.frame( + time = c(4, 5), + event = c(1, 1), + y = c(4, 3), + z = c(2, 4), + x = c(NA, NA), + w = c("h", "i") + ) + suppressMessages({ + actual <- remove_missing_rows(df, Surv(time, event) ~ y, extra_vars = c("z", "w")) + }) + + rownames(expected) <- NULL + rownames(actual) <- NULL + + expect_equal(expected, actual) }) test_that("expand_initial_values smoke tests", { - vals <- list("a" = 1, "b" = 2, "c" = c(1, 2)) - siz <- list("a" = 5, "b" = 1, "c" = 2) - actual <- expand_initial_values(vals, siz) - expected <- list("a" = c(1, 1, 1, 1, 1), "b" = 2, "c" = c(1, 2)) - expect_equal(actual, expected) + vals <- list("a" = 1, "b" = 2, "c" = c(1, 2)) + siz <- list("a" = 5, "b" = 1, "c" = 2) + actual <- expand_initial_values(vals, siz) + expected <- list("a" = c(1, 1, 1, 1, 1), "b" = 2, "c" = c(1, 2)) + expect_equal(actual, expected) }) test_that("replace_with_lookup smoke tests", { - - vals <- list(1, "b", "a", 4, 5) - lku <- list("a" = 3, "b" = 2, "c" = 4) - actual <- replace_with_lookup(vals, lku) - expected <- list(1, 2, 3, 4, 5) - expect_equal(actual, expected) - - - vals <- list(1, "b", "a", 4, c(5, 6)) - lku <- list("a" = 3, "b" = 2, "c" = 4) - expect_error( - replace_with_lookup(vals, lku), - regexp = "`sizes` must be length 1 numerics" - ) - + vals <- list(1, "b", "a", 4, 5) + lku <- list("a" = 3, "b" = 2, "c" = 4) + actual <- replace_with_lookup(vals, lku) + expected <- list(1, 2, 3, 4, 5) + expect_equal(actual, expected) + + + vals <- list(1, "b", "a", 4, c(5, 6)) + lku <- list("a" = 3, "b" = 2, "c" = 4) + expect_error( + replace_with_lookup(vals, lku), + regexp = "`sizes` must be length 1 numerics" + ) }) - - - - From 5b0595d60de1f632a3be292b33e9b4e7a1568fe9 Mon Sep 17 00:00:00 2001 From: cicdguy <26552821+cicdguy@users.noreply.github.com> Date: Wed, 10 May 2023 09:10:41 -0500 Subject: [PATCH 8/8] Run workflows