inst/doc/longitudinal.R

In pcvr: Plant Phenotyping and Bayesian Statistics

## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)

## ----message = FALSE----------------------------------------------------------
library(pcvr)
library(data.table) # for fread
library(ggplot2)
library(patchwork) # for easy ggplot manipulation/combination
library(brms)

## ----eval = FALSE-------------------------------------------------------------
# if (!"cmdstanr" %in% installed.packages()) {
#   install.packages("cmdstanr", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))
# }
# if (!"brms" %in% installed.packages()) {
#   install.packages("brms")
# }
# library(brms)
# library(cmdstanr)
# cmdstanr::install_cmdstan()

## -----------------------------------------------------------------------------
simdf <- growthSim("logistic", n = 20, t = 25, params = list(
  "A" = c(200, 160),
  "B" = c(13, 11),
  "C" = c(3, 3.5)
))
l <- ggplot(simdf, aes(time, y, group = interaction(group, id))) +
  geom_line(aes(color = group)) +
  labs(title = "Logistic") +
  theme_minimal() +
  theme(legend.position = "none")

simdf <- growthSim("gompertz", n = 20, t = 25, params = list(
  "A" = c(200, 160),
  "B" = c(13, 11),
  "C" = c(0.2, 0.25)
))
g <- ggplot(simdf, aes(time, y, group = interaction(group, id))) +
  geom_line(aes(color = group)) +
  labs(title = "Gompertz") +
  theme_minimal() +
  theme(legend.position = "none")

simdf <- growthSim("monomolecular", n = 20, t = 25, params = list("A" = c(200, 160), "B" = c(0.08, 0.1)))
m <- ggplot(simdf, aes(time, y, group = interaction(group, id))) +
  geom_line(aes(color = group)) +
  labs(title = "Monomolecular") +
  theme_minimal() +
  theme(legend.position = "none")

simdf <- growthSim("exponential", n = 20, t = 25, params = list("A" = c(15, 20), "B" = c(0.095, 0.095)))
e <- ggplot(simdf, aes(time, y, group = interaction(group, id))) +
  geom_line(aes(color = group)) +
  labs(title = "Exponential") +
  theme_minimal() +
  theme(legend.position = "none")

simdf <- growthSim("linear", n = 20, t = 25, params = list("A" = c(1.1, 0.95)))
ln <- ggplot(simdf, aes(time, y, group = interaction(group, id))) +
  geom_line(aes(color = group)) +
  labs(title = "Linear") +
  theme_minimal() +
  theme(legend.position = "none")

simdf <- growthSim("power law", n = 20, t = 25, params = list("A" = c(16, 11), "B" = c(0.75, 0.7)))
pl <- ggplot(simdf, aes(time, y, group = interaction(group, id))) +
  geom_line(aes(color = group)) +
  labs(title = "Power Law") +
  theme_minimal() +
  theme(legend.position = "none")

patch <- (l + g + m) / (e + ln + pl)
patch

## -----------------------------------------------------------------------------
set.seed(345)
gomp <- growthSim("gompertz", n = 20, t = 35, params = list(
  "A" = c(200, 180, 160),
  "B" = c(20, 22, 18),
  "C" = c(0.15, 0.2, 0.1)
))


sigma_df <- aggregate(y ~ group + time, data = gomp, FUN = sd)

ggplot(sigma_df, aes(x = time, y = y, group = group)) +
  geom_line(color = "gray60") +
  pcv_theme() +
  labs(y = "SD of y", title = "Gompertz Sigma")

## ----message = FALSE----------------------------------------------------------
draw_gomp_sigma <- function(x) {
  return(23 * exp(-21 * exp(-0.22 * x)))
}
ggplot(sigma_df, aes(x = time, y = y)) +
  geom_line(aes(group = group), color = "gray60") +
  geom_hline(aes(yintercept = 12, color = "Homoskedastic"),
    linetype = 5,
    key_glyph = draw_key_path
  ) +
  geom_abline(aes(slope = 0.8, intercept = 0, color = "Linear"),
    linetype = 5,
    key_glyph = draw_key_path
  ) +
  geom_smooth(
    method = "gam", aes(color = "Spline"), linetype = 5, se = FALSE,
    key_glyph = draw_key_path
  ) +
  geom_function(fun = draw_gomp_sigma, aes(color = "Gompertz"), linetype = 5) +
  scale_color_viridis_d(option = "plasma", begin = 0.1, end = 0.9) +
  guides(color = guide_legend(override.aes = list(linewidth = 1, linetype = 1))) +
  pcv_theme() +
  theme(legend.position = "bottom") +
  labs(y = "SD of y", title = "Gompertz Sigma", color = "")

## -----------------------------------------------------------------------------
ss <- growthSS(
  model = "gompertz", form = y ~ time | id / group, sigma = "int",
  df = gomp, start = list("A" = 130, "B" = 15, "C" = 0.25)
)

ss

## ----eval=FALSE---------------------------------------------------------------
# fit_h <- fitGrowth(ss, iter = 1000, cores = 4, chains = 4, silent = 0)
# 
# brmPlot(fit_h, form = ss$pcvrForm, df = ss$df)

## -----------------------------------------------------------------------------
ss <- growthSS(
  model = "gompertz", form = y ~ time | id / group, sigma = "linear",
  df = gomp, start = list("A" = 130, "B" = 15, "C" = 0.25)
)

ss

## ----eval=FALSE---------------------------------------------------------------
# fit_l <- fitGrowth(ss,
#   iter = 1000, cores = 4, chains = 4, silent = 0,
#   control = list(adapt_delta = 0.999, max_treedepth = 20)
# )
# 
# p1 <- brmPlot(fit_l, form = ss$pcvrForm, df = ss$df)
# p2 <- p1 + coord_cartesian(ylim = c(0, 300))
# p <- p1 / p2
# p

## -----------------------------------------------------------------------------
ss <- growthSS(
  model = "gompertz", form = y ~ time | id / group, sigma = "spline",
  df = gomp, start = list("A" = 130, "B" = 15, "C" = 0.25)
)

ss

## ----eval=FALSE---------------------------------------------------------------
# fit_s <- fitGrowth(ss,
#   iter = 2000, cores = 4, chains = 4, silent = 0,
#   control = list(adapt_delta = 0.999, max_treedepth = 20)
# )
# 
# brmPlot(fit_s, form = ss$pcvrForm, df = ss$df)

## -----------------------------------------------------------------------------
ss <- growthSS(
  model = "gompertz", form = y ~ time | id / group, sigma = "gompertz",
  df = gomp, start = list(
    "A" = 130, "B" = 15, "C" = 0.25,
    "sigmaA" = 15, "sigmaB" = 15, "sigmaC" = 0.25
  ),
  type = "brms"
)

ss

## ----eval=FALSE---------------------------------------------------------------
# fit_g <- fitGrowth(ss,
#   iter = 2000, cores = 4, chains = 4, silent = 0,
#   control = list(adapt_delta = 0.999, max_treedepth = 20)
# )
# 
# brmPlot(fit_g, form = ss$pcvrForm, df = ss$df)

## ----message = FALSE----------------------------------------------------------
draw_gomp_sigma <- function(x) {
  return(23 * exp(-21 * exp(-0.22 * x)))
}
draw_logistic_sigma <- function(x) {
  return(20 / (1 + exp((15 - x) / 2)))
}
draw_logistic_exp <- function(x) {
  return(2.5 * exp(0.08 * x))
}
draw_logistic_quad <- function(x) {
  return((0.3 * x) + (0.02 * x^2))
}

ggplot(sigma_df, aes(x = time, y = y)) +
  geom_line(aes(group = group), color = "gray60", linetype = 5) +
  geom_hline(aes(yintercept = 12, color = "Homoskedastic"), linetype = 1) +
  geom_abline(aes(slope = 0.8, intercept = 0, color = "Linear"),
    linetype = 1,
    key_glyph = draw_key_path
  ) +
  geom_smooth(
    method = "gam", aes(color = "Spline"), linetype = 1, se = FALSE,
    key_glyph = draw_key_path
  ) +
  geom_function(fun = draw_gomp_sigma, aes(color = "Gompertz"), linetype = 1) +
  geom_function(fun = draw_logistic_sigma, aes(color = "Logistic"), linetype = 1) +
  geom_function(fun = draw_logistic_exp, aes(color = "Exponential"), linetype = 1) +
  geom_function(fun = draw_logistic_quad, aes(color = "Quadratic"), linetype = 1) +
  scale_color_viridis_d(option = "plasma", begin = 0.1, end = 0.9) +
  guides(color = guide_legend(override.aes = list(linewidth = 1, linetype = 1))) +
  pcv_theme() +
  theme(legend.position = "bottom") +
  labs(y = "SD of y", title = "Gompertz Sigma", color = "")

## ----eval = FALSE-------------------------------------------------------------
# loo_spline <- add_criterion(fit_s, "loo")
# loo_homo <- add_criterion(fit_h, "loo")
# loo_linear <- add_criterion(fit_l, "loo")
# loo_gomp <- add_criterion(fit_g, "loo")
# 
# h <- loo_homo$criteria$loo$estimates[3, 1]
# s <- loo_spline$criteria$loo$estimates[3, 1]
# l <- loo_linear$criteria$loo$estimates[3, 1]
# g <- loo_gomp$criteria$loo$estimates[3, 1]
# 
# loodf <- data.frame(loo = c(h, s, l, g), model = c("Homosked", "Spline", "Linear", "Gompertz"))
# loodf$model <- factor(loodf$model, levels = unique(loodf$model[order(loodf$loo)]), ordered = TRUE)
# 
# ggplot(
#   loodf,
#   aes(x = model, y = loo, fill = model)
# ) +
#   geom_col() +
#   scale_fill_viridis_d() +
#   labs(y = "LOO Information Criteria", x = "Sub Model of Sigma") +
#   theme_minimal() +
#   theme(legend.position = "none")

## ----eval = FALSE-------------------------------------------------------------
# set.seed(345)
# ln <- growthSim("linear", n = 5, t = 10, params = list("A" = c(2, 3, 10)))
# 
# strongPrior <- prior(student_t(3, 0, 5), dpar = "sigma", class = "b") +
#   prior(gamma(2, 0.1), class = "nu", lb = 0.001) +
#   prior(normal(10, .05), nlpar = "A", lb = 0)
# 
# ss <- growthSS(
#   model = "linear", form = y ~ time | id / group, sigma = "homo",
#   df = ln, priors = strongPrior
# )
# 
# fit <- fitGrowth(ss, iter = 1000, cores = 2, chains = 2, silent = 0)
# 
# brmPlot(fit, form = ss$pcvrForm, df = ss$df) +
#   coord_cartesian(ylim = c(0, 100))

## ----eval = FALSE-------------------------------------------------------------
# weakPrior <- prior(student_t(3, 0, 5), dpar = "sigma", class = "b") +
#   prior(gamma(2, 0.1), class = "nu", lb = 0.001) +
#   prior(lognormal(log(10), 0.25), nlpar = "A", lb = 0)
# 
# ss <- growthSS(
#   model = "linear", form = y ~ time | id / group, sigma = "homo",
#   df = ln, priors = weakPrior
# )
# 
# fit <- fitGrowth(ss, iter = 1000, cores = 2, chains = 2, silent = 0)
# 
# brmPlot(fit, form = ss$pcvrForm, df = ss$df) +
#   coord_cartesian(ylim = c(0, 100))

## -----------------------------------------------------------------------------
priors <- list("A" = 130, "B" = 10, "C" = 0.2)
priorPlots <- plotPrior(priors)
priorPlots[[1]] / priorPlots[[2]] / priorPlots[[3]]

## -----------------------------------------------------------------------------
twoPriors <- list("A" = c(100, 130), "B" = c(6, 12), "C" = c(0.5, 0.25))
plotPrior(twoPriors, "gompertz", n = 100)[[1]]

## -----------------------------------------------------------------------------
set.seed(123)
mv_df <- mvSim(dists = list(rnorm = list(mean = 100, sd = 30)), wide = FALSE)
mv_df$group <- rep(c("a", "b"), times = 900)
mv_df <- mv_df[mv_df$value > 0, ]
mv_df$label <- as.numeric(gsub("sim_", "", mv_df$variable))

ss1 <- mvSS(
  model = "linear", form = label | value ~ group, df = mv_df,
  start = list("A" = 5), type = "brms", spectral_index = "ci_rededge"
)

ss1

## -----------------------------------------------------------------------------
set.seed(123)
simdf <- growthSim("logistic", n = 20, t = 25, params = list(
  "A" = c(200, 160),
  "B" = c(13, 11),
  "C" = c(3, 3.5)
))

## -----------------------------------------------------------------------------
nls_ss <- growthSS(
  model = "logistic", form = y ~ time | id / group,
  df = simdf, type = "nls"
)

## -----------------------------------------------------------------------------
nlrq_ss <- growthSS(
  model = "logistic", form = y ~ time | id / group,
  df = simdf, type = "nlrq",
  tau = seq(0.01, 0.99, 0.04)
)

## -----------------------------------------------------------------------------
nlme_ss <- growthSS(
  model = "logistic", form = y ~ time | id / group,
  df = simdf, sigma = "power", type = "nlme"
)

## -----------------------------------------------------------------------------
mgcv_ss <- growthSS(
  model = "gam", form = y ~ time | id / group,
  df = simdf, type = "mgcv"
)

## ----eval=FALSE---------------------------------------------------------------
# brms_ss <- growthSS(
#   model = "logistic", form = y ~ time | id / group,
#   sigma = "spline", df = simdf,
#   start = list("A" = 130, "B" = 10, "C" = 1)
# )

## -----------------------------------------------------------------------------
ggplot(simdf, aes(time, y, group = interaction(group, id))) +
  geom_line(aes(color = group))

## ----warning=FALSE, message=FALSE---------------------------------------------
nls_fit <- fitGrowth(nls_ss)
nlrq_fit <- fitGrowth(nlrq_ss)
nlme_fit <- fitGrowth(nlme_ss)
mgcv_fit <- fitGrowth(mgcv_ss)

## ----eval=FALSE---------------------------------------------------------------
# brms_fit <- fitGrowth(brms_ss,
#   iter = 500, cores = 1, chains = 1,
#   control = list(adapt_delta = 0.999, max_treedepth = 20)
# )

## ----warning=FALSE, message=FALSE---------------------------------------------
growthPlot(nls_fit, form = nls_ss$pcvrForm, df = nls_ss$df)
growthPlot(nlrq_fit, form = nlrq_ss$pcvrForm, df = nlrq_ss$df)
growthPlot(nlme_fit, form = nlme_ss$pcvrForm, df = nlme_ss$df)
growthPlot(mgcv_fit, form = mgcv_ss$pcvrForm, df = mgcv_ss$df)

## ----eval=FALSE---------------------------------------------------------------
# growthPlot(brms_fit, form = brms_ss$pcvrForm, df = brms_ss$df)

## -----------------------------------------------------------------------------
testGrowth(nls_ss, nls_fit, test = "A")$anova

## -----------------------------------------------------------------------------
testGrowth(nlrq_ss, nlrq_fit, test = "A")[["0.49"]]

## ----warning=FALSE, message=FALSE---------------------------------------------
testGrowth(nlme_ss, nlme_fit, test = "A")$anova

## -----------------------------------------------------------------------------
testGrowth(mgcv_ss, mgcv_fit)$anova

## -----------------------------------------------------------------------------
testGrowth(fit = nls_fit, test = list(
  "A1 - A2 *1.1",
  "(B1+1) - B2",
  "C1 - (C2-0.5)",
  "A1/B1 - (1.1 * A2/B2)"
))

## -----------------------------------------------------------------------------
testGrowth(fit = nlme_fit, test = list(
  "(A.groupa / A.groupb) - 0.9",
  "1 + (B.groupa - B.groupb)",
  "C.groupa/C.groupb - 1"
))

## ----echo=TRUE, eval=FALSE----------------------------------------------------
# (hyp <- brms::hypothesis(brms_fit, "(A_groupa) > 1.1 * (A_groupb)"))

## ----echo=FALSE, eval=TRUE----------------------------------------------------
structure(list(
  Hypothesis = "((A_groupa))-(1.1*(A_groupb)) > 0",
  Estimate = 16.7880224, Est.Error = 3.13518501598807, CI.Lower = 12.012705,
  CI.Upper = 21.952505, Evid.Ratio = Inf, Post.Prob = 1, Star = "*"
), row.names = c(NA, -1L), class = "data.frame")

## ----echo=TRUE, eval=FALSE----------------------------------------------------
# simdf <- growthSim(
#   model = "linear + linear",
#   n = 20, t = 25,
#   params = list("linear1A" = c(15, 12), "changePoint1" = c(8, 6), "linear2A" = c(3, 5))
# )
# 
# ss <- growthSS(
#   model = "linear + linear", form = y ~ time | id / group, sigma = "spline",
#   start = list("linear1A" = 10, "changePoint1" = 5, "linear2A" = 2),
#   df = simdf, type = "brms"
# )
# 
# fit <- fitGrowth(ss, backend = "cmdstanr", iter = 500, chains = 1, cores = 1)

## ----echo=TRUE, eval=FALSE----------------------------------------------------
# simdf <- growthSim("linear + logistic",
#   n = 20, t = 25,
#   params = list(
#     "linear1A" = c(15, 12), "changePoint1" = c(8, 6),
#     "logistic2A" = c(100, 150), "logistic2B" = c(10, 8),
#     "logistic2C" = c(3, 2.5)
#   )
# )
# 
# ss <- growthSS(
#   model = "linear + logistic", form = y ~ time | id / group, sigma = "spline",
#   list(
#     "linear1A" = 10, "changePoint1" = 5,
#     "logistic2A" = 100, "logistic2B" = 10, "logistic2C" = 3
#   ),
#   df = simdf, type = "brms"
# )
# 
# fit <- fitGrowth(ss, backend = "cmdstanr", iter = 500, chains = 1, cores = 1)

## ----eval=FALSE, echo=TRUE----------------------------------------------------
# ss <- growthSS(
#   model = "linear + gam", form = y ~ time | id / group, sigma = "int",
#   list("linear1A" = 10, "changePoint1" = 5),
#   df = simdf, type = "brms"
# )
# 
# fit <- fitGrowth(ss, backend = "cmdstanr", iter = 500, chains = 1, cores = 1)

## ----echo=TRUE, eval=FALSE----------------------------------------------------
# simdf <- growthSim("linear + linear + linear",
#   n = 25, t = 50,
#   params = list(
#     "linear1A" = c(10, 12), "changePoint1" = c(8, 6),
#     "linear2A" = c(1, 2), "changePoint2" = c(25, 30), "linear3A" = c(20, 24)
#   )
# )
# 
# ss <- growthSS(
#   model = "linear + linear + linear", form = y ~ time | id / group, sigma = "spline",
#   list(
#     "linear1A" = 10, "changePoint1" = 5,
#     "linear2A" = 2, "changePoint2" = 15,
#     "linear3A" = 5
#   ), df = simdf, type = "brms"
# )
# 
# fit <- fitGrowth(ss, backend = "cmdstanr", iter = 500, chains = 1, cores = 1)

## ----echo=FALSE, eval=FALSE---------------------------------------------------
# set.seed(123)
# noise <- do.call(rbind, lapply(1:30, function(i) {
#   chngpt <- rnorm(2, 18, 2)
#   noise_i <- rbind(
#     data.frame(
#       id = paste0("id_", i), time = 1:chngpt[1], group = "a",
#       y = c(runif(chngpt[1] - 1, 0, 20), rnorm(1, 5, 1))
#     ),
#     data.frame(
#       id = paste0("id_", i), time = 1:chngpt[2], group = "b",
#       y = c(runif(chngpt[2] - 1, 0, 20), rnorm(1, 5, 1))
#     )
#   )
#   return(noise_i)
# }))
# noise2 <- do.call(rbind, lapply(1:30, function(i) {
#   start1 <- max(noise[noise$id == paste0("id_", i) & noise$group == "a", "time"])
#   start2 <- max(noise[noise$id == paste0("id_", i) & noise$group == "b", "time"])
#   noise2_i <- rbind(
#     data.frame(
#       id = paste0("id_", i), time = start1:40, group = "a",
#       y = c(runif(length(start1:40), 15, 50))
#     ),
#     data.frame(
#       id = paste0("id_", i), time = start2:40, group = "b",
#       y = c(runif(length(start2:40), 15, 50))
#     )
#   )
#   return(noise2_i)
# }))
# simdf <- rbind(noise, noise2)

## ----echo=TRUE, eval=FALSE----------------------------------------------------
# ss <- growthSS(
#   model = "int + int", form = y ~ time | id / group, sigma = "int + int",
#   list(
#     "int1" = 10, "changePoint1" = 10, "int2" = 20, # main model
#     "sigmaint1" = 10, "sigmachangePoint1" = 10, "sigmaint2" = 10
#   ), # sub model
#   df = simdf, type = "brms"
# )
# 
# fit <- fitGrowth(ss, backend = "cmdstanr", iter = 500, chains = 1, cores = 1)

## ----echo=FALSE, eval=FALSE---------------------------------------------------
# set.seed(123)
# noise <- do.call(rbind, lapply(1:30, function(i) {
#   chngpt <- rnorm(2, 18, 2)
#   noise_i <- rbind(
#     data.frame(
#       id = paste0("id_", i), time = 1:chngpt[1], group = "a",
#       y = c(runif(chngpt[1] - 1, 0, 20), rnorm(1, 5, 1))
#     ),
#     data.frame(
#       id = paste0("id_", i), time = 1:chngpt[2], group = "b",
#       y = c(runif(chngpt[2] - 1, 0, 20), rnorm(1, 5, 1))
#     )
#   )
#   return(noise_i)
# }))
# signal <- growthSim("linear",
#   n = 30, t = 20,
#   params = list("A" = c(3, 5))
# )
# signal <- do.call(rbind, lapply(unique(paste0(signal$id, signal$group)), function(int) {
#   noisesub <- noise[paste0(noise$id, noise$group) == int, ]
#   signalSub <- signal[paste0(signal$id, signal$group) == int, ]
#   y_end <- noisesub[noisesub$time == max(noisesub$time), "y"]
#   signalSub$time <- signalSub$time + max(noisesub$time)
#   signalSub$y <- y_end + signalSub$y
#   return(signalSub)
# }))
# simdf <- rbind(noise, signal)

## ----echo=TRUE, eval=FALSE----------------------------------------------------
# ss <- growthSS(
#   model = "int + linear", form = y ~ time | id / group, sigma = "int + linear",
#   list(
#     "int1" = 10, "changePoint1" = 10, "linear2A" = 20,
#     "sigmaint1" = 10, "sigmachangePoint1" = 10, "sigmalinear2A" = 10
#   ),
#   df = simdf, type = "brms"
# )
# fit <- fitGrowth(ss, backend = "cmdstanr", iter = 500, chains = 1, cores = 1)

## ----eval=FALSE, echo=FALSE---------------------------------------------------
# set.seed(123)
# noise <- do.call(rbind, lapply(1:30, function(i) {
#   chngpt <- rnorm(2, 18, 2)
#   noise_i <- rbind(
#     data.frame(
#       id = paste0("id_", i), time = 1:chngpt[1], group = "a",
#       y = c(runif(chngpt[1] - 1, 0, 20), rnorm(1, 5, 1))
#     ),
#     data.frame(
#       id = paste0("id_", i), time = 1:chngpt[2], group = "b",
#       y = c(runif(chngpt[2] - 1, 0, 20), rnorm(1, 5, 1))
#     )
#   )
#   return(noise_i)
# }))
# signal <- growthSim("logistic",
#   n = 20, t = 30,
#   params = list("A" = c(200, 160), "B" = c(13, 11), "C" = c(3, 3.5))
# )
# signal <- do.call(rbind, lapply(unique(paste0(signal$id, signal$group)), function(int) {
#   noisesub <- noise[paste0(noise$id, noise$group) == int, ]
#   signalSub <- signal[paste0(signal$id, signal$group) == int, ]
#   y_end <- noisesub[noisesub$time == max(noisesub$time), "y"]
#   signalSub$time <- signalSub$time + max(noisesub$time)
#   signalSub$y <- y_end + signalSub$y
#   return(signalSub)
# }))
# simdf <- rbind(noise, signal)
# simdf <- simdf[simdf$time < 45, ]

## ----echo=TRUE, eval=FALSE----------------------------------------------------
# ss <- growthSS(
#   model = "int+logistic", form = y ~ time | id / group, sigma = "int + spline",
#   list(
#     "int1" = 5, "changePoint1" = 10,
#     "logistic2A" = 130, "logistic2B" = 10, "logistic2C" = 3,
#     "sigmaint1" = 5, "sigmachangePoint1" = 15
#   ),
#   df = simdf, type = "brms"
# )
# fit <- fitGrowth(ss, backend = "cmdstanr", iter = 500, chains = 1, cores = 1)

## ----eval=FALSE---------------------------------------------------------------
# print(load(url("https://raw.githubusercontent.com/joshqsumner/pcvrTestData/main/brmsFits.rdata")))
# from3to25 <- list(
#   fit_3, fit_5, fit_7, fit_9, fit_11, fit_13,
#   fit_15, fit_17, fit_19, fit_21, fit_23, fit_25
# )
# 
# distributionPlot(fits = from3to25, form = y ~ time | id / group, params = c("A", "B", "C"), d = simdf)