diff --git a/.Rbuildignore b/.Rbuildignore index af495e10..ba6e45da 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -14,4 +14,5 @@ ^pkgdown$ ^vignettes/articles$ ^DO_NOT_COMMIT* -^\.#.*$ \ No newline at end of file +^demon* +^#.*#$ diff --git a/.gitignore b/.gitignore index c01668a4..8bf3c4f6 100644 --- a/.gitignore +++ b/.gitignore @@ -15,3 +15,4 @@ docs tools/lib tools/*.html \#\*R:scripts\*\# +*.org diff --git a/CRAN-SUBMISSION b/CRAN-SUBMISSION index ba784e12..1052c881 100644 --- a/CRAN-SUBMISSION +++ b/CRAN-SUBMISSION @@ -1,3 +1,3 @@ -Version: 1.3.1 -Date: 2025-08-18 20:38:57 UTC -SHA: 03f0cc66b17667d39907be289ad0930cfb6ea1f6 +Version: 1.4.1 +Date: 2026-06-15 15:59:18 UTC +SHA: dbea2a109738231b67030e826fee70e27ca53408 diff --git a/DESCRIPTION b/DESCRIPTION index 1a139044..879f50cd 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: pcvr Type: Package Title: Plant Phenotyping and Bayesian Statistics -Version: 1.4.0 +Version: 1.4.1 Authors@R: c(person("Josh", "Sumner", email = "jsumner@danforthcenter.org", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-3399-9063")), @@ -21,8 +21,7 @@ License: GPL-2 Encoding: UTF-8 LazyData: true LazyDataCompression: xz -Additional_repositories: https://mc-stan.org/r-packages/ -RoxygenNote: 7.3.2 +Additional_repositories: https://stan-dev.r-universe.dev Imports: FactoMineR, rlang, @@ -64,3 +63,5 @@ VignetteBuilder: knitr Config/testthat/edition: 3 URL: https://github.com/danforthcenter/pcvr, https://plantcv.org/, https://danforthcenter.github.io/pcvr/ BugReports: https://github.com/danforthcenter/pcvr/issues +Config/roxygen2/version: 8.0.0 +RoxygenNote: 7.3.2 diff --git a/NEWS.md b/NEWS.md index b82cb9ba..8d6f38d4 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,7 @@ +# pcvr 1.4.1 + +Minor edits in tutorials/articles/examples for CRAN and accessibility. + # pcvr 1.4.0 Added `stat_growthss` to make ggplot layers of models from `growthSS`. diff --git a/R/barg.R b/R/barg.R index a40d766a..bcc3c8da 100644 --- a/R/barg.R +++ b/R/barg.R @@ -103,25 +103,6 @@ #' @seealso \link{plotPrior} for visual prior predictive checks. #' @examples #' \donttest{ -#' simdf <- growthSim("logistic", -#' n = 20, t = 25, -#' params = list("A" = c(200, 160), "B" = c(13, 11), "C" = c(3, 3.5)) -#' ) -#' ss <- growthSS( -#' model = "logistic", form = y ~ time | id / group, sigma = "logistic", -#' df = simdf, start = list( -#' "A" = 130, "B" = 12, "C" = 3, -#' "sigmaA" = 20, "sigmaB" = 10, "sigmaC" = 2 -#' ), type = "brms" -#' ) -#' fit_test <- fitGrowth(ss, -#' iter = 600, cores = 1, chains = 1, backend = "cmdstanr", -#' sample_prior = "only" # only sampling from prior for speed -#' ) -#' b <- barg(fit_test, ss) -#' fit_2 <- fit_test -#' fit_list <- list(fit_test, fit_2) -#' x <- barg(fit_list, list(ss, ss)) #' #' x <- conjugate( #' s1 = rnorm(10, 10, 1), s2 = rnorm(10, 13, 1.5), method = "t", diff --git a/R/brmPlot.R b/R/brmPlot.R index 377d5d40..2ca048b8 100644 --- a/R/brmPlot.R +++ b/R/brmPlot.R @@ -26,17 +26,7 @@ #' @importFrom stats as.formula #' @examples #' \donttest{ -#' simdf <- growthSim( -#' "logistic", -#' n = 20, t = 25, -#' params = list("A" = c(200, 160), "B" = c(13, 11), "C" = c(3, 3.5)) -#' ) -#' ss <- growthSS( -#' model = "logistic", form = y ~ time | id / group, sigma = "spline", -#' list("A" = 130, "B" = 10, "C" = 3), -#' df = simdf, type = "brms" -#' ) -#' fit <- fitGrowth(ss, backend = "cmdstanr", iter = 500, chains = 1, cores = 1) +#' data(fit) #' growthPlot(fit = fit, form = y ~ time | group, groups = "a", df = ss$df) #' } #' diff --git a/R/brmSurvPlot.R b/R/brmSurvPlot.R index 4643061c..bd204c77 100644 --- a/R/brmSurvPlot.R +++ b/R/brmSurvPlot.R @@ -24,21 +24,12 @@ #' n = 20, t = 50, #' params = list("A" = c(1, 1), "B" = c(0.15, 0.1)) #' ) -#' ss1 <- growthSS( +#' ss <- growthSS( #' model = "survival weibull", form = y > 100 ~ time | id / group, #' df = df, start = c(0, 5) #' ) -#' fit1 <- fitGrowth(ss1, iter = 600, cores = 2, chains = 2, backend = "cmdstanr") -#' brmSurvPlot(fit1, form = ss1$pcvrForm, df = ss1$df) -#' -#' # note that using the cumulative hazard to calculate survival is likely to underestimate -#' # survival in these plots if events do not start immediately. -#' ss2 <- growthSS( -#' model = "survival binomial", form = y > 100 ~ time | id / group, -#' df = df, start = c(-4, 3) -#' ) -#' fit2 <- fitGrowth(ss2, iter = 600, cores = 2, chains = 2, backend = "cmdstanr") -#' brmSurvPlot(fit2, form = ss2$pcvrForm, df = ss2$df) +#' data(surv) +#' brmSurvPlot(surv, form = ss$pcvrForm, df = ss$df) #' } #' #' @return Returns a ggplot showing a brms model's credible diff --git a/R/brmViolin.R b/R/brmViolin.R index 744a1dd4..5bee7b3e 100644 --- a/R/brmViolin.R +++ b/R/brmViolin.R @@ -21,18 +21,7 @@ #' @examples #' \donttest{ #' set.seed(123) -#' simdf <- growthSim( -#' "logistic", -#' n = 20, t = 25, -#' params = list("A" = c(200, 180, 190, 160), "B" = c(13, 11, 10, 10), "C" = c(3, 3, 3.25, 3.5)) -#' ) -#' ss <- growthSS( -#' model = "logistic", form = y ~ time | id / group, sigma = "int", -#' list("A" = 130, "B" = 10, "C" = 3), -#' df = simdf, type = "brms" -#' ) -#' -#' fit <- fitGrowth(ss, backend = "cmdstanr", iter = 500, chains = 1, cores = 1) +#' data(fit) #' brmViolin(fit, ss, ".../A_groupd > 1.05") # all groups used #' brmViolin(fit, ss, "abs(1 - ((...)/(C_groupd - B_groupd))) > 0.05") # rather arbitrary #' brmViolin(fit, ss, "abs(1 - ((...)/(C_groupa - B_groupd))) > 0.05") # totally arbitrary diff --git a/R/combineDraws.R b/R/combineDraws.R index 4c111e15..cba905b0 100644 --- a/R/combineDraws.R +++ b/R/combineDraws.R @@ -24,61 +24,12 @@ #' @examples #' # note that this example will fit several models using Stan and may run slowly. #' \donttest{ -#' simdf <- growthSim("logistic", -#' n = 20, t = 25, -#' params = list( -#' "A" = c(200, 160, 220, 200, 140, 300), -#' "B" = c(13, 11, 10, 9, 16, 12), -#' "C" = c(3, 3.5, 3.2, 2.8, 3.3, 2.5) -#' ) -#' ) -#' ss_ab <- growthSS( -#' model = "logistic", form = y ~ time | id / group, -#' sigma = "logistic", df = simdf[simdf$group %in% c("a", "b"), ], -#' start = list( -#' "A" = 130, "B" = 12, "C" = 3, -#' "sigmaA" = 15, "sigmaB" = 10, "sigmaC" = 3 -#' ), type = "brms" -#' ) -#' -#' ss_cd <- growthSS( -#' model = "logistic", form = y ~ time | id / group, -#' sigma = "logistic", df = simdf[simdf$group %in% c("c", "d"), ], -#' start = list( -#' "A" = 130, "B" = 12, "C" = 3, -#' "sigmaA" = 15, "sigmaB" = 10, "sigmaC" = 3 -#' ), type = "brms" -#' ) -#' -#' ss_ef <- growthSS( -#' model = "logistic", form = y ~ time | id / group, -#' sigma = "logistic", df = simdf[simdf$group %in% c("e", "f"), ], -#' start = list( -#' "A" = 130, "B" = 12, "C" = 3, -#' "sigmaA" = 15, "sigmaB" = 10, "sigmaC" = 3 -#' ), type = "brms" -#' ) -#' ss_ef2 <- growthSS( -#' model = "gompertz", form = y ~ time | id / group, -#' sigma = "logistic", df = simdf[simdf$group %in% c("e", "f"), ], -#' start = list( -#' "A" = 130, "B" = 12, "C" = 3, -#' "sigmaA" = 15, "sigmaB" = 10, "sigmaC" = 3 -#' ), type = "brms" -#' ) -#' -#' -#' fit_ab <- fitGrowth(ss_ab, chains = 1, cores = 1, iter = 1000) -#' fit_ab2 <- fitGrowth(ss_ab, chains = 1, cores = 1, iter = 1200) -#' fit_cd <- fitGrowth(ss_cd, chains = 1, cores = 1, iter = 1000) -#' fit_ef <- fitGrowth(ss_ef, chains = 1, cores = 1, iter = 1000) -#' fit_ef2 <- fitGrowth(ss_ef2, chains = 1, cores = 1, iter = 1000) -#' -#' x <- combineDraws(fit_ab, fit_cd, fit_ef) -#' draws_ef <- as.data.frame(fit_ef) -#' draws_ef <- draws_ef[, grepl("^b_", colnames(draws_ef))] -#' x2 <- combineDraws(fit_ab2, fit_cd, draws_ef) -#' x3 <- combineDraws(fit_ab, fit_cd, fit_ef2) +#' data(fit) +#' fit_1 <- fit_2 <- fit +#' x <- combineDraws(fit_1, fit_2, fit) +#' draws_1 <- as.data.frame(fit_1) +#' draws_1 <- draws_1[, grepl("^b_", colnames(draws_ef))] +#' x2 <- combineDraws(fit_2, draws_1) #' } #' #' @return Returns a dataframe of posterior draws. diff --git a/R/datasets.R b/R/datasets.R new file mode 100644 index 00000000..f4daebfe --- /dev/null +++ b/R/datasets.R @@ -0,0 +1,26 @@ +#' Example model +#' +#' @description A logistic brms model of simulated growth data. +#' +#' @format A brmsfit object and a growthss object +#' @examples +#' \donttest{ +#' data(fit, package = "pcvr") +#' summary(fit) +#' } +#' +"fit" + + +#' Example survival model +#' +#' @description An exponential survival brms model of simulated data. +#' +#' @format A brmsfit object and a growthss object +#' @examples +#' \donttest{ +#' data(surv, package = "pcvr") +#' summary(surv) +#' } +#' +"surv" diff --git a/R/growthSS.R b/R/growthSS.R index 65eaaa57..74ec5519 100644 --- a/R/growthSS.R +++ b/R/growthSS.R @@ -336,15 +336,7 @@ #' ) #' lapply(ss, class) #' ss$initfun() -#' # the next step would typically be compiling/fitting the model -#' # here we use very few chains and very few iterations for speed, but more of both is better. -#' \donttest{ -#' fit_test <- fitGrowth(ss, -#' iter = 500, cores = 1, chains = 1, backend = "cmdstanr", -#' control = list(adapt_delta = 0.999, max_treedepth = 20) -#' ) -#' } -#' +#' # the next step would typically be compiling/fitting the model with fitGrowth #' #' # formulas and priors will look different if there is only one group in the data #' diff --git a/R/mvSS.R b/R/mvSS.R index add75639..fed2bf4b 100644 --- a/R/mvSS.R +++ b/R/mvSS.R @@ -34,12 +34,6 @@ #' model = "linear", form = label | value ~ group, df = mv_df, #' start = list("A" = 5), type = "brms", spectral_index = "none" #' ) -#' \donttest{ -#' mod1 <- fitGrowth(ss1, backend = "cmdstanr", iter = 1000, chains = 1, cores = 1) -#' growthPlot(mod1, ss1$pcvrForm, df = ss1$df) -#' library(ggplot2) -#' ggplot() + stat_brms_model(fit = mod1, ss = ss1) -#' } #' #' # when the model is longitudinal the same model is possible with growthSS #' @@ -97,13 +91,6 @@ #' } #' })) #' -#' \donttest{ -#' if (rlang::is_installed("mnormt")) { -#' m2 <- fitGrowth(ss_mv1, backend = "cmdstanr", iter = 1000, chains = 1, cores = 1) -#' growthPlot(m2, ss_mv1$pcvrForm, df = ss_mv1$df) -#' } -#' } -#' #' @export mvSS <- function(model = "linear", form, sigma = NULL, df, start = NULL, diff --git a/R/stat_growthSS.R b/R/stat_growthSS.R index 673855cc..2245be33 100644 --- a/R/stat_growthSS.R +++ b/R/stat_growthSS.R @@ -33,15 +33,7 @@ #' #' @examples #' library(ggplot2) -#' simdf <- growthSim("logistic", -#' n = 20, t = 25, -#' params = list("A" = c(200, 160), "B" = c(13, 11), "C" = c(3, 3.5)) -#' ) -#' ss <- growthSS( -#' model = "logistic", form = y ~ time | id / group, -#' df = simdf, start = NULL, type = "nls" -#' ) -#' fit <- fitGrowth(ss) +#' data(fit) #' ggplot() + #' stat_growthss(fit = fit, ss = ss) #' diff --git a/cran-comments.md b/cran-comments.md index 5dcd6dc1..7b38c09d 100644 --- a/cran-comments.md +++ b/cran-comments.md @@ -1,42 +1,36 @@ -# pcvr 1.3.1 +# pcvr 1.4.1 ## R CMD check results -0 errors ✔ | 0 warnings ✔ | 2 notes ✖ +0 errors ✔ | 0 warnings ✔ | 1 notes ✖ -❯ checking CRAN incoming feasibility ... [4s/6s] NOTE - Maintainer: ‘Josh Sumner ’ - - Possibly misspelled words in DESCRIPTION: - Kruschke (16:2, 17:2, 17:46) - Phenotyping (3:14) - phenotyping (12:44) +❯ checking for future file timestamps ... NOTE + unable to verify current time +## Notes + +Resubmitting due to cran testing error probably related to either not installing +suggested packages or the change from `mc-stan.org` to `stan-dev.r--universe.dev` +for some additional packages and for development features +(ggplot2 geom support, more distributions, etc). +Previously a NOTE related to the failure specified: + +``` +❯ checking package dependencies ... NOTE +Package suggested but not available for checking: ‘cmdstanr’ Suggests or Enhances not in mainstream repositories: cmdstanr Availability using Additional_repositories specification: - cmdstanr yes https://mc-stan.org/r-packages/ -❯ checking package namespace information ... OK -❯ checking package dependencies ... NOTE -Package suggested but not available for checking: ‘cmdstanr’ - -## Notes + cmdstanr yes https://stan-dev.r-universe.dev +``` -Resubmitting due to large changes in the development version. +But the `cmdstanr` suggestion is available from the new additional repository specified in DESCRIPTION. No notes seem critical. + Names (PlantCV, Kruschke) and words (phenotyping) in DESCRIPTION are not misspelled. -The `cmdstanr` suggestion is available from the additional repository specified in DESCRIPTION. The R CMD check on MacOs via github actions yields a NOTE about installed package size as well. -``` -* checking installed package size ... NOTE - installed size is 6.4Mb - sub-directories of 1Mb or more: - doc 4.1Mb - R 2.0Mb -``` - Additionally there are some examples wrapped in `\donttest` that may take several minutes to run. ## Test environments diff --git a/data/fit.rda b/data/fit.rda new file mode 100644 index 00000000..247d73e8 Binary files /dev/null and b/data/fit.rda differ diff --git a/data/surv.rda b/data/surv.rda new file mode 100644 index 00000000..d9f602f4 Binary files /dev/null and b/data/surv.rda differ diff --git a/man/barg.Rd b/man/barg.Rd index ad453617..9af27aef 100644 --- a/man/barg.Rd +++ b/man/barg.Rd @@ -110,25 +110,6 @@ by the provided list. If you are unsure what random-generation function to use t } \examples{ \donttest{ -simdf <- growthSim("logistic", - n = 20, t = 25, - params = list("A" = c(200, 160), "B" = c(13, 11), "C" = c(3, 3.5)) -) -ss <- growthSS( - model = "logistic", form = y ~ time | id / group, sigma = "logistic", - df = simdf, start = list( - "A" = 130, "B" = 12, "C" = 3, - "sigmaA" = 20, "sigmaB" = 10, "sigmaC" = 2 - ), type = "brms" -) -fit_test <- fitGrowth(ss, - iter = 600, cores = 1, chains = 1, backend = "cmdstanr", - sample_prior = "only" # only sampling from prior for speed -) -b <- barg(fit_test, ss) -fit_2 <- fit_test -fit_list <- list(fit_test, fit_2) -x <- barg(fit_list, list(ss, ss)) x <- conjugate( s1 = rnorm(10, 10, 1), s2 = rnorm(10, 13, 1.5), method = "t", diff --git a/man/brmPlot.Rd b/man/brmPlot.Rd index 918083c2..04129eab 100644 --- a/man/brmPlot.Rd +++ b/man/brmPlot.Rd @@ -51,17 +51,7 @@ means) can be visualized easily using this function. This will generally be call } \examples{ \donttest{ -simdf <- growthSim( - "logistic", - n = 20, t = 25, - params = list("A" = c(200, 160), "B" = c(13, 11), "C" = c(3, 3.5)) -) -ss <- growthSS( - model = "logistic", form = y ~ time | id / group, sigma = "spline", - list("A" = 130, "B" = 10, "C" = 3), - df = simdf, type = "brms" -) -fit <- fitGrowth(ss, backend = "cmdstanr", iter = 500, chains = 1, cores = 1) +data(fit) growthPlot(fit = fit, form = y ~ time | group, groups = "a", df = ss$df) } diff --git a/man/brmSurvPlot.Rd b/man/brmSurvPlot.Rd index c4a436ca..2e8b6826 100644 --- a/man/brmSurvPlot.Rd +++ b/man/brmSurvPlot.Rd @@ -46,21 +46,12 @@ df <- growthSim("exponential", n = 20, t = 50, params = list("A" = c(1, 1), "B" = c(0.15, 0.1)) ) -ss1 <- growthSS( +ss <- growthSS( model = "survival weibull", form = y > 100 ~ time | id / group, df = df, start = c(0, 5) ) -fit1 <- fitGrowth(ss1, iter = 600, cores = 2, chains = 2, backend = "cmdstanr") -brmSurvPlot(fit1, form = ss1$pcvrForm, df = ss1$df) - -# note that using the cumulative hazard to calculate survival is likely to underestimate -# survival in these plots if events do not start immediately. -ss2 <- growthSS( - model = "survival binomial", form = y > 100 ~ time | id / group, - df = df, start = c(-4, 3) -) -fit2 <- fitGrowth(ss2, iter = 600, cores = 2, chains = 2, backend = "cmdstanr") -brmSurvPlot(fit2, form = ss2$pcvrForm, df = ss2$df) +data(surv) +brmSurvPlot(surv, form = ss$pcvrForm, df = ss$df) } } diff --git a/man/brmViolin.Rd b/man/brmViolin.Rd index 68f20446..392a8d62 100644 --- a/man/brmViolin.Rd +++ b/man/brmViolin.Rd @@ -31,18 +31,7 @@ Function to visualize hypotheses tested on brms models similar to those made usi \examples{ \donttest{ set.seed(123) -simdf <- growthSim( - "logistic", - n = 20, t = 25, - params = list("A" = c(200, 180, 190, 160), "B" = c(13, 11, 10, 10), "C" = c(3, 3, 3.25, 3.5)) -) -ss <- growthSS( - model = "logistic", form = y ~ time | id / group, sigma = "int", - list("A" = 130, "B" = 10, "C" = 3), - df = simdf, type = "brms" -) - -fit <- fitGrowth(ss, backend = "cmdstanr", iter = 500, chains = 1, cores = 1) +data(fit) brmViolin(fit, ss, ".../A_groupd > 1.05") # all groups used brmViolin(fit, ss, "abs(1 - ((...)/(C_groupd - B_groupd))) > 0.05") # rather arbitrary brmViolin(fit, ss, "abs(1 - ((...)/(C_groupa - B_groupd))) > 0.05") # totally arbitrary diff --git a/man/combineDraws.Rd b/man/combineDraws.Rd index 144efcba..3464dbc3 100644 --- a/man/combineDraws.Rd +++ b/man/combineDraws.Rd @@ -37,61 +37,12 @@ the list by default so passing the \code{names} argument is helpful and can be d \examples{ # note that this example will fit several models using Stan and may run slowly. \donttest{ -simdf <- growthSim("logistic", - n = 20, t = 25, - params = list( - "A" = c(200, 160, 220, 200, 140, 300), - "B" = c(13, 11, 10, 9, 16, 12), - "C" = c(3, 3.5, 3.2, 2.8, 3.3, 2.5) - ) -) -ss_ab <- growthSS( - model = "logistic", form = y ~ time | id / group, - sigma = "logistic", df = simdf[simdf$group \%in\% c("a", "b"), ], - start = list( - "A" = 130, "B" = 12, "C" = 3, - "sigmaA" = 15, "sigmaB" = 10, "sigmaC" = 3 - ), type = "brms" -) - -ss_cd <- growthSS( - model = "logistic", form = y ~ time | id / group, - sigma = "logistic", df = simdf[simdf$group \%in\% c("c", "d"), ], - start = list( - "A" = 130, "B" = 12, "C" = 3, - "sigmaA" = 15, "sigmaB" = 10, "sigmaC" = 3 - ), type = "brms" -) - -ss_ef <- growthSS( - model = "logistic", form = y ~ time | id / group, - sigma = "logistic", df = simdf[simdf$group \%in\% c("e", "f"), ], - start = list( - "A" = 130, "B" = 12, "C" = 3, - "sigmaA" = 15, "sigmaB" = 10, "sigmaC" = 3 - ), type = "brms" -) -ss_ef2 <- growthSS( - model = "gompertz", form = y ~ time | id / group, - sigma = "logistic", df = simdf[simdf$group \%in\% c("e", "f"), ], - start = list( - "A" = 130, "B" = 12, "C" = 3, - "sigmaA" = 15, "sigmaB" = 10, "sigmaC" = 3 - ), type = "brms" -) - - -fit_ab <- fitGrowth(ss_ab, chains = 1, cores = 1, iter = 1000) -fit_ab2 <- fitGrowth(ss_ab, chains = 1, cores = 1, iter = 1200) -fit_cd <- fitGrowth(ss_cd, chains = 1, cores = 1, iter = 1000) -fit_ef <- fitGrowth(ss_ef, chains = 1, cores = 1, iter = 1000) -fit_ef2 <- fitGrowth(ss_ef2, chains = 1, cores = 1, iter = 1000) - -x <- combineDraws(fit_ab, fit_cd, fit_ef) -draws_ef <- as.data.frame(fit_ef) -draws_ef <- draws_ef[, grepl("^b_", colnames(draws_ef))] -x2 <- combineDraws(fit_ab2, fit_cd, draws_ef) -x3 <- combineDraws(fit_ab, fit_cd, fit_ef2) +data(fit) +fit_1 <- fit_2 <- fit +x <- combineDraws(fit_1, fit_2, fit) +draws_1 <- as.data.frame(fit_1) +draws_1 <- draws_1[, grepl("^b_", colnames(draws_ef))] +x2 <- combineDraws(fit_2, draws_1) } } diff --git a/man/fit.Rd b/man/fit.Rd new file mode 100644 index 00000000..0d06f4a6 --- /dev/null +++ b/man/fit.Rd @@ -0,0 +1,23 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/datasets.R +\docType{data} +\name{fit} +\alias{fit} +\title{Example model} +\format{ +A brmsfit object and a growthss object +} +\usage{ +fit +} +\description{ +A logistic brms model of simulated growth data. +} +\examples{ +\donttest{ +data(fit, package = "pcvr") +summary(fit) +} + +} +\keyword{datasets} diff --git a/man/growthSS.Rd b/man/growthSS.Rd index 6b84647e..4761b8d4 100644 --- a/man/growthSS.Rd +++ b/man/growthSS.Rd @@ -356,15 +356,7 @@ ss <- growthSS( ) lapply(ss, class) ss$initfun() -# the next step would typically be compiling/fitting the model -# here we use very few chains and very few iterations for speed, but more of both is better. -\donttest{ -fit_test <- fitGrowth(ss, - iter = 500, cores = 1, chains = 1, backend = "cmdstanr", - control = list(adapt_delta = 0.999, max_treedepth = 20) -) -} - +# the next step would typically be compiling/fitting the model with fitGrowth # formulas and priors will look different if there is only one group in the data diff --git a/man/mvSS.Rd b/man/mvSS.Rd index f775f9fe..ba1b6f74 100644 --- a/man/mvSS.Rd +++ b/man/mvSS.Rd @@ -67,12 +67,6 @@ ss1 <- mvSS( model = "linear", form = label | value ~ group, df = mv_df, start = list("A" = 5), type = "brms", spectral_index = "none" ) -\donttest{ -mod1 <- fitGrowth(ss1, backend = "cmdstanr", iter = 1000, chains = 1, cores = 1) -growthPlot(mod1, ss1$pcvrForm, df = ss1$df) -library(ggplot2) -ggplot() + stat_brms_model(fit = mod1, ss = ss1) -} # when the model is longitudinal the same model is possible with growthSS @@ -130,13 +124,6 @@ unlist(lapply(names(ss_mv1), function(nm) { } })) -\donttest{ -if (rlang::is_installed("mnormt")) { - m2 <- fitGrowth(ss_mv1, backend = "cmdstanr", iter = 1000, chains = 1, cores = 1) - growthPlot(m2, ss_mv1$pcvrForm, df = ss_mv1$df) -} -} - } \seealso{ \link{fitGrowth} for fitting the model specified by this list. diff --git a/man/stat_growthss.Rd b/man/stat_growthss.Rd index 83fb70ec..200f8341 100644 --- a/man/stat_growthss.Rd +++ b/man/stat_growthss.Rd @@ -107,15 +107,7 @@ The geometries used for each type of model are: } \examples{ library(ggplot2) -simdf <- growthSim("logistic", - n = 20, t = 25, - params = list("A" = c(200, 160), "B" = c(13, 11), "C" = c(3, 3.5)) -) -ss <- growthSS( - model = "logistic", form = y ~ time | id / group, - df = simdf, start = NULL, type = "nls" -) -fit <- fitGrowth(ss) +data(fit) ggplot() + stat_growthss(fit = fit, ss = ss) diff --git a/man/surv.Rd b/man/surv.Rd new file mode 100644 index 00000000..f80d20fb --- /dev/null +++ b/man/surv.Rd @@ -0,0 +1,23 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/datasets.R +\docType{data} +\name{surv} +\alias{surv} +\title{Example survival model} +\format{ +A brmsfit object and a growthss object +} +\usage{ +surv +} +\description{ +An exponential survival brms model of simulated data. +} +\examples{ +\donttest{ +data(surv, package = "pcvr") +summary(surv) +} + +} +\keyword{datasets} diff --git a/tools/make_data.R b/tools/make_data.R new file mode 100644 index 00000000..88c9ba73 --- /dev/null +++ b/tools/make_data.R @@ -0,0 +1,24 @@ +setwd("~/pcvr/data/") +devtools::load_all("~/pcvr") +set.seed(123) +df <- pcvr::growthSim("logistic", n = 10, t = 20, params = list(A = 100, B = 10, C = 3)) +ss <- pcvr::growthSS("logistic", y ~ time | id / group, sigma = "logistic", df = df, + start = list(A = 100, B = 10, C = 3, sigmaA = 10, sigmaB = 10, sigmaC = 3)) +fit <- pcvr::fitGrowth(ss, backend = "cmdstanr", iter = 4100, warmup = 4000, chains = 2, cores = 2) +names(fit$fit) +object.size(fit) #165128 bytes +object.size(ss) # 35368 bytes +save(fit, ss, file = "fit.rda") +# might want a survival one? +set.seed(123) +df <- growthSim("exponential", + n = 10, t = 30, + params = list("A" = c(1, 1), "B" = c(0.15, 0.2)) +) +ss <- growthSS( + model = "survival weibull", form = y > 25 ~ time | id / group, + df = df, start = c(0, 5) +) +fit <- fitGrowth(ss, iter = 4100, warmup = 4000, cores = 2, chains = 2, backend = "cmdstanr") +brmSurvPlot(fit, form = ss$pcvrForm, df = ss$df) +save(fit, ss, file = "surv.rda") diff --git a/vignettes/articles/field_capacity.Rmd b/vignettes/articles/field_capacity.Rmd index 72e49f6a..d610d26e 100644 --- a/vignettes/articles/field_capacity.Rmd +++ b/vignettes/articles/field_capacity.Rmd @@ -584,6 +584,8 @@ head(d) Target field capacities are in the data that the phenotyping core provides and observed field capacities can be calculated from the watering metadata. ```{r} +#| fig.alt: > +#| Example of field capacity in a drought treatment that is working well set.seed(123) d2 <- data.frame( day = rep(1:20, 4), @@ -619,7 +621,11 @@ ggplot(d[d$variable == "field_capacity", ], ) + pcv_theme() + theme(legend.position = "none") +``` +```{r} +#| fig.alt: > +#| Example of field capacity in a drought treatment that is not working well ggplot(d[d$variable == "field_capacity", ], aes(x = day, y = value)) + facet_wrap(~treatment) + diff --git a/vignettes/articles/reading_pcv_data.Rmd b/vignettes/articles/reading_pcv_data.Rmd index 0c8d8b84..0f348c77 100644 --- a/vignettes/articles/reading_pcv_data.Rmd +++ b/vignettes/articles/reading_pcv_data.Rmd @@ -163,7 +163,6 @@ out <- pcv.time(sv, plantingDelay = 0, phenotype = "area_pixels", cutoff = 10, timeCol = "timestamp", group = c("barcode", "rotation"), plot = TRUE ) -out$plot sv <- out$data dim(sv) ``` @@ -173,6 +172,8 @@ Note that in these plots, and particularly clearly in the DAE plot, there are lo Before moving on we'll also check the grouping in our data. Here we see that we have lots of plants with more than one image per day. ```{r} +#| fig.alt: > +#| Checking grouping of data checkGroups(sv, c("DAS", "barcode", "rotation", "genotype", "fertilizer")) ``` @@ -194,6 +195,8 @@ dim(sv_ag_with_outliers) The `pcv.outliers` function can be used to remove outliers relative to a phenotype using cook's distance. Here due to the experimental design having these plants germinate on the machine around 1 percent of data are removed as outliers. The plot shows removed data points in red, although here that is hard to see. ```{r} +#| fig.alt: > +#| Plotting outliers out <- pcv.outliers( df = sv_ag_with_outliers, phenotype = "area_pixels", group = c("DAS", "genotype", "fertilizer"), plotgroup = c("barcode") @@ -206,6 +209,8 @@ dim(sv_ag) It is also useful to check our grouping assumptions again, here we see that there are some plants with multiple images from a single day. ```{r} +#| fig.alt: > +#| Checking Groups again after aggregation/outlier removal checkGroups(sv_ag, c("DAS", "barcode", "genotype", "fertilizer")) ``` @@ -215,6 +220,8 @@ checkGroups(sv_ag, c("DAS", "barcode", "genotype", "fertilizer")) We might also want to check the watering data, which can be read easily from json with `pcv.water`. ```{r} +#| fig.alt: > +#| Plotting watering data water <- pcv.water(paste0(base_url, "metadata.json")) water$genotype <- substr(water$barcode, 3, 5) @@ -244,6 +251,8 @@ A common use for watering data is to look at water use efficiency (WUE). Here we ```{r} +#| fig.alt: > +#| Plotting Pseudo water use efficiency test <- pwue(df = sv_ag, w = water, pheno = "area_pixels", time = "timestamp", id = "barcode") ggplot(test, aes(x = DAS, y = pWUE, color = genotype, group = barcode)) +