Skip to content

Commit 192f217

Browse files
committed
Merge branch 'main' of https://github.com/UCD-SERG/serodynamics into Kwan_2025
# Conflicts: # .Rbuildignore # .gitignore # DESCRIPTION # NEWS.md # tests/testthat/_snaps/run_mod/nostrat-curve-params.csv # tests/testthat/_snaps/run_mod/sim-strat-curve-params.csv # tests/testthat/_snaps/run_mod/strat-curve-params.csv
2 parents b327aa3 + 4ee4b99 commit 192f217

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

41 files changed

+1535024
-240
lines changed

.Rbuildignore

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -22,4 +22,7 @@
2222
^\.project$
2323
^\.settings$
2424
^\.lintr\.R$
25+
^_extensions$
26+
^vignettes/\.quarto$
27+
^vignettes/*_files$
2528

.gitignore

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -21,4 +21,5 @@ serodynamics*.tgz
2121

2222
inst/doc
2323
docs
24+
**/.quarto/
2425

.lintr.R

Lines changed: 21 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11

22

3-
undesirable_functions <-
4-
lintr::default_undesirable_functions |>
3+
undesirable_functions <-
4+
lintr::default_undesirable_functions |>
55
lintr::modify_defaults(
66

77
# following https://github.com/r-lib/devtools/blob/2aa51ef/.lintr.R:
@@ -25,25 +25,26 @@ undesirable_functions <-
2525
"\nSee:\n",
2626
"<https://r-pkgs.org/code.html#sec-code-r-landscape> and\n",
2727
"<https://r-pkgs.org/testing-design.html#sec-testing-design-self-contained>",
28-
"\nfor more details"
28+
"\nfor more details."
2929
),
3030

31-
structure = NULL
32-
# see https://github.com/r-lib/lintr/pull/2227 and
31+
structure = NULL,
32+
browser = NULL
33+
# see https://github.com/r-lib/lintr/pull/2227 and
3334
# rebuttal https://github.com/r-lib/lintr/pull/2227#issuecomment-1800302675
3435

3536
)
3637

3738
# define snake_case with uppercase acronyms allowed;
3839
# see https://github.com/r-lib/lintr/issues/2844 for details:
3940
withr::local_package("rex")
40-
snake_case_ACRO = rex::rex(
41+
snake_case_ACROs1 <- rex::rex(
4142
start,
4243
maybe("."),
43-
some_of(lower, digit) %or% some_of(upper, digit),
44+
list(some_of(upper), maybe("s"), zero_or_more(digit)) %or% list(some_of(lower), zero_or_more(digit)),
4445
zero_or_more(
4546
"_",
46-
some_of(lower, digit) %or% some_of(upper, digit)
47+
list(some_of(upper), maybe("s"), zero_or_more(digit)) %or% list(some_of(lower), zero_or_more(digit))
4748
),
4849
end
4950
)
@@ -54,7 +55,7 @@ linters <- lintr::linters_with_defaults(
5455
lintr::redundant_equals_linter(),
5556
lintr::pipe_consistency_linter(pipe = "|>"),
5657
lintr::object_name_linter(
57-
regexes = c(snake_case_ACRO = snake_case_ACRO)
58+
regexes = c(snake_case_ACROs1 = snake_case_ACROs1)
5859
),
5960
lintr::undesirable_function_linter(
6061
fun = undesirable_functions,
@@ -64,9 +65,19 @@ linters <- lintr::linters_with_defaults(
6465

6566
# prevent warnings from lintr::read_settings:
6667
rm(undesirable_functions)
67-
rm(snake_case_ACRO)
68+
rm(snake_case_ACROs1)
6869
exclusions <- list(
6970
`data-raw` = list(
7071
pipe_consistency_linter = Inf
72+
),
73+
vignettes = list(
74+
undesirable_function_linter = Inf
75+
),
76+
"inst/analyses" = list(
77+
undesirable_function_linter = Inf
78+
),
79+
"tests/testthat.R" = list(
80+
undesirable_function_linter = Inf
7181
)
82+
7283
)

DESCRIPTION

Lines changed: 4 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,7 @@ Authors@R: c(
55
person("Peter", "Teunis", , "p.teunis@emory.edu", role = c("aut", "cph"),
66
comment = "Author of the method and original code."),
77
person("Samuel", "Schildhauer", , "sschildhauer@ucdavis.edu", role = c("aut", "cre")),
8-
person("Kwan Ho", "Lee", role = "ctb"),
8+
person("Kwan Ho", "Lee", , "ksjlee@ucdavis.edu", role = "ctb"),
99
person("Kristen", "Aiemjoy", , "kaiemjoy@ucdavis.edu", role = "aut"),
1010
person("Douglas Ezra", "Morrison", , "demorrison@ucdavis.edu", role = "aut")
1111
)
@@ -48,7 +48,8 @@ Suggests:
4848
rmarkdown,
4949
magrittr,
5050
ssdtools,
51-
rex
51+
rex,
52+
quarto
5253
Language: en-US
5354
Config/testthat/edition: 3
5455
Config/needs/test:
@@ -67,4 +68,4 @@ Remotes:
6768
Depends:
6869
R (>= 4.1.0)
6970
LazyData: true
70-
VignetteBuilder: knitr
71+
VignetteBuilder: knitr, quarto

NEWS.md

Lines changed: 8 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -4,8 +4,15 @@
44

55
## New features
66

7-
87
* Added `plot_predicted_curve()` (#68)
8+
* Replacing old data object with new run_mod output (#102)
9+
* Adding class assignment to run_mod output (#76)
10+
* Making prep_priors modifiable (#78)
11+
* Changes to `run_mod()` output:
12+
- Taking out `include_subs` as an input option, default will include all
13+
individuals
14+
- Making a single tbl as output
15+
- All other pieces will be attributes.
916
* Changes to `run_mod()` (#79):
1017
- `jags.post` now optionally included in output, as specified by argument
1118
`with_post`

R/Run_Mod.R

Lines changed: 60 additions & 50 deletions
Original file line numberDiff line numberDiff line change
@@ -1,20 +1,20 @@
11
#' @title Run Jags Model
22
#' @author Sam Schildhauer
33
#' @description
4-
#' run_mod() takes a data frame and adjustable mcmc inputs to
5-
#' [runjags::run.jags()] as an mcmc
6-
#' bayesian model to estimate antibody dynamic curve parameters.
4+
#' `run_mod()` takes a data frame and adjustable MCMC inputs to
5+
#' [runjags::run.jags()] as an MCMC
6+
#' Bayesian model to estimate antibody dynamic curve parameters.
77
#' The [rjags::jags.model()] models seroresponse dynamics to an
88
#' infection. The antibody dynamic curve includes the following parameters:
99
#' - y0 = baseline antibody concentration
1010
#' - y1 = peak antibody concentration
1111
#' - t1 = time to peak
12-
#' - r = shape parameter
12+
#' - shape = shape parameter
1313
#' - alpha = decay rate
1414
#' @param data A [base::data.frame()] with the following columns.
1515
#' @param file_mod The name of the file that contains model structure.
1616
#' @param nchain An [integer] between 1 and 4 that specifies
17-
#' the number of mcmc chains to be run per jags model.
17+
#' the number of MCMC chains to be run per jags model.
1818
#' @param nadapt An [integer] specifying the number of adaptations per chain.
1919
#' @param nburn An [integer] specifying the number of burn ins before sampling.
2020
#' @param nmc An [integer] specifying the number of samples in posterior chains.
@@ -26,43 +26,50 @@
2626
#' should be included as an element of the [list] object returned by `run_mod()`
2727
#' (see `Value` section below for details).
2828
#' Note: These objects can be large.
29-
#' @param include_subs A [logical] value specifying whether posterior
30-
#' distributions should be included for all subjects. A value of [FALSE] will
31-
#' only include the predictive distribution.
32-
#' @returns A [list] containing the following elements:
33-
#' - `"jags.post"`: a [list] containing one or more [runjags::runjags-class]
34-
#' objects (one per stratum).
35-
#' - A [base::data.frame()] titled `curve_params` that contains the posterior
36-
#' distribution will be exported with the following attributes:
37-
#' - `iteration` = number of sampling iterations
38-
#' - `chain` = number of mcmc chains run; between 1 and 4
39-
#' - `indexid` = "newperson", indicating posterior distribution
40-
#' - `antigen_iso` = antibody/antigen type combination being evaluated
41-
#' - `alpha` = posterior estimate of decay rate
42-
#' - `r` = posterior estimate of shape parameter
43-
#' - `t1` = posterior estimate of time to peak
44-
#' - `y0` = posterior estimate of baseline antibody concentration
45-
#' - `y1` = posterior estimate of peak antibody concentration
46-
#' - `stratified variable` = the variable used to stratify jags model
47-
#' - A [list] of `attributes` that summarize the jags inputs, including:
29+
#' @returns An `sr_model` class object: a subclass of [dplyr::tbl_df] that
30+
#' contains MCMC samples from the joint posterior distribution of the model
31+
#' parameters, conditional on the provided input `data`,
32+
#' including the following:
33+
#' - `iteration` = Number of sampling iterations
34+
#' - `chain` = Number of MCMC chains run; between 1 and 4
35+
#' - `Parameter` = Parameter being estimated. Includes the following:
36+
#' - `y0` = Posterior estimate of baseline antibody concentration
37+
#' - `y1` = Posterior estimate of peak antibody concentration
38+
#' - `t1` = Posterior estimate of time to peak
39+
#' - `shape` = Posterior estimate of shape parameter
40+
#' - `alpha` = Posterior estimate of decay rate
41+
#' - `Iso_type` = Antibody/antigen type combination being evaluated
42+
#' - `Stratification` = The variable used to stratify jags model
43+
#' - `Subject` = ID of subject being evaluated
44+
#' - `value` = Estimated value of the parameter
45+
#' - The following [attributes] are included in the output:
4846
#' - `class`: Class of the output object.
4947
#' - `nChain`: Number of chains run.
5048
#' - `nParameters`: The amount of parameters estimated in the model.
5149
#' - `nIterations`: Number of iteration specified.
5250
#' - `nBurnin`: Number of burn ins.
53-
#' - `nThin`: Thinning number (niter/nmc).
51+
#' - `nThin`: Thinning number (niter/nmc)
52+
#' - `priors`: A [list] that summarizes the input priors, including:
53+
#' - `mu_hyp_param`
54+
#' - `prec_hyp_param`
55+
#' - `omega_param`
56+
#' - `wishdf`
57+
#' - `prec_logy_hyp_param`
58+
#' - An optional `"jags.post"` attribute, included when argument
59+
#' `with_post` = TRUE.
60+
#' @inheritDotParams prep_priors
5461
#' @export
5562
#' @example inst/examples/run_mod-examples.R
5663
run_mod <- function(data,
57-
file_mod,
64+
file_mod = serodynamics_example("model.jags"),
5865
nchain = 4,
5966
nadapt = 0,
6067
nburn = 0,
6168
nmc = 100,
6269
niter = 100,
6370
strat = NA,
6471
with_post = FALSE,
65-
include_subs = FALSE) {
72+
...) {
6673
## Conditionally creating a stratification list to loop through
6774
if (is.na(strat)) {
6875
strat_list <- "None"
@@ -97,7 +104,8 @@ run_mod <- function(data,
97104

98105
# prepare data for modeline
99106
longdata <- prep_data(dl_sub)
100-
priors <- prep_priors(max_antigens = longdata$n_antigen_isos)
107+
priorspec <- prep_priors(max_antigens = longdata$n_antigen_isos,
108+
...)
101109

102110
# inputs for jags model
103111
nchains <- nchain # nr of MC chains to run simultaneously
@@ -111,7 +119,7 @@ run_mod <- function(data,
111119

112120
jags_post <- runjags::run.jags(
113121
model = file_mod,
114-
data = c(longdata, priors),
122+
data = c(longdata, priorspec),
115123
inits = initsfunction,
116124
method = "parallel",
117125
adapt = nadapt,
@@ -127,13 +135,13 @@ run_mod <- function(data,
127135
# stratification and will only be included if specified.
128136
jags_post_final[[i]] <- jags_post
129137

130-
# Unpacking and cleaning mcmc output.
138+
# Unpacking and cleaning MCMC output.
131139
jags_unpack <- ggmcmc::ggs(jags_post[["mcmc"]])
132140

133141
# Adding attributes
134142
mod_atts <- attributes(jags_unpack)
135-
# Only keeping necesarry attributes
136-
mod_atts <- mod_atts[3:8]
143+
# Only keeping necessary attributes
144+
mod_atts <- mod_atts[4:8]
137145

138146
# extracting antigen-iso combinations to correctly number
139147
# then by the order they are estimated by the program.
@@ -169,25 +177,27 @@ run_mod <- function(data,
169177
jags_out <- jags_out[complete.cases(jags_out), ]
170178
# Outputting the finalized jags output as a data frame with the
171179
# jags output results for each stratification rbinded.
172-
# Logical argument to include posterior of all subjects or just the
173-
# predictive distribution (new person).
174-
if (!include_subs) {
175-
jags_out <- jags_out |>
176-
filter(.data$Subject == "newperson")
177-
}
178180

179-
# Logical argument to keep the raw jags post or not.
181+
# Making output a tibble and restructing.
182+
jags_out <- dplyr::as_tibble(jags_out) |>
183+
select(!c("Parameter")) |>
184+
rename("Parameter" = "Parameter_sub")
185+
jags_out <- jags_out[, c("Iteration", "Chain", "Parameter", "Iso_type",
186+
"Stratification", "Subject", "value")]
187+
current_atts <- attributes(jags_out)
188+
current_atts <- c(current_atts, mod_atts)
189+
attributes(jags_out) <- current_atts
190+
191+
# Adding priors
192+
jags_out <- jags_out |>
193+
structure("priors" = attributes(priorspec)$used_priors)
194+
195+
# Conditionally adding jags.post
180196
if (with_post) {
181-
jags_out <- list(
182-
"curve_params" = jags_out,
183-
"jags.post" = jags_post_final,
184-
"attributes" = mod_atts
185-
)
186-
} else {
187-
jags_out <- list(
188-
"curve_params" = jags_out,
189-
"attributes" = mod_atts
190-
)
191-
}
197+
jags_out <- jags_out |>
198+
structure(jags.post = jags_post_final)
199+
}
200+
jags_out <- jags_out |>
201+
structure(class = union("sr_model", class(jags_out)))
192202
jags_out
193203
}

R/nepal_sees_jags_output.R

Lines changed: 16 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -8,15 +8,25 @@
88
#' which is the diagnosis type (typhoid or
99
#' paratyphoid). Keeping only IDs `"newperson"`, `"sees_npl_1"`, `"sees_npl_2"`.
1010
#'
11-
#' @format ## `nepal_sees_jags_output`
12-
#' A [list] consisting of the following named elements:
13-
#' \describe{
14-
#' \item{curve_params}{A [data.frame] titled `curve_params` that contains the
11+
#' @format An S3 object of class `sr_model`: A [dplyr::tbl_df] that contains the
1512
#' posterior predictive distribution of the person-specific parameters for a
1613
#' "new person" with no observed data (`Subject = "newperson"`) and posterior
1714
#' distributions of the person-specific parameters for two arbitrarily-chosen
18-
#' subjects (`"sees_npl_1"` and `"sees_npl_2"`)}
19-
#' \item{attributes}{A [list] of `attributes` that summarize the jags inputs}
15+
#' subjects (`"sees_npl_1"` and `"sees_npl_2"`).
16+
#' Contains 40,000 `rows`, 7 `columns`, and model `attributes`.
17+
#' \describe{
18+
#' \item{Iteration}{Number of sampling iterations: 500 iterations}
19+
#' \item{Chain}{Number of MCMC chains run: 2 chains run}
20+
#' \item{Parameter}{Parameter being estimated}
21+
#' \item{Iso_type}{Antibody/antigen type combination being evaluated:
22+
#' `HlyE_IgA` and `HlyE_IgG`}
23+
#' \item{Stratification}{The variable used to stratify jags model: `typhi` and
24+
#' `paratyphi`}
25+
#' \item{Subject}{ID of subject being evaluated: `newperson`, `sees_npl_1`,
26+
#' `sees_npl_2`}
27+
#' \item{value}{Estimated value of the parameter}
28+
#' \item{attributes}{A [list] of `attributes` that summarize the jags inputs,
29+
#' priors, and optional jags_post mcmc object}
2030
#' }
2131
#' @source reference study: <https://doi.org/10.1016/S2666-5247(22)00114-8>
2232
"nepal_sees_jags_output"

R/plot_jags_densitydx.R

Lines changed: 7 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -32,16 +32,16 @@
3232
#' @example inst/examples/examples-plot_jags_densitydx.R
3333

3434
plot_jags_dens <- function(data,
35-
iso = unique(data$curve_params$Iso_type),
36-
param = unique(data$curve_params$Parameter_sub),
37-
strat = unique(data$curve_params$Stratification)) {
38-
visualize_jags <- data[["curve_params"]]
35+
iso = unique(data$Iso_type),
36+
param = unique(data$Parameter),
37+
strat = unique(data$Stratification)) {
38+
3939
attributes_jags <- data[["attributes"]]
4040

4141
dens_strat_list <- list()
4242
for (i in strat) {
4343

44-
visualize_jags_sub <- visualize_jags |>
44+
visualize_jags_sub <- data |>
4545
dplyr::filter(.data$Stratification == i) |>
4646
dplyr::filter(.data$Subject == "newperson")
4747

@@ -55,12 +55,12 @@ plot_jags_dens <- function(data,
5555
# Will not loop through parameters, as we may want each to show on the
5656
# same plot by default.
5757
visualize_jags_plot <- visualize_jags_plot |>
58-
dplyr::filter(.data$Parameter_sub %in% param)
58+
dplyr::filter(.data$Parameter %in% param)
5959

6060
visualize_jags_plot <- visualize_jags_plot |>
6161
# Changing parameter name to reflect the input
6262
dplyr::mutate(Parameter = paste0("iso = ", j, ", parameter = ",
63-
.data$Parameter_sub, ", strat = ",
63+
.data$Parameter, ", strat = ",
6464
i),
6565
value = log(.data$value))
6666
# Assigning attributes, which are needed to run ggs_density

0 commit comments

Comments
 (0)