Package 'pcvr'

Description

subset helper function for use reading in large data, called in pcv.sub.read

Usage

awkHelper(inputFile, filters, awk = NULL)

Arguments

Details

awkHelper attempts to make awk commands from human readable input. Currently when filters are supplied the input file has quotes removed by 'sed' then is piped into awk, so an equivalent command line statement may be: sed 's/\"//g' pcvrTest2.csv | awk -F ',' '{ if (NR==1 || $18=="area") { print } }'

Value

Returns a character string representing a unix style awk statement which is typically passed to pipe or used as a connection in data.table::fread.

Examples

tryCatch(
  { # in case offline
    link1 <- "https://gist.githubusercontent.com/seankross/"
    link2 <- "a412dfbd88b3db70b74b/raw/5f23f993cd87c283ce766e7ac6b329ee7cc2e1d1/mtcars.csv"
    file <- paste0(link1, link2)
    awkHelper(file, list("gear in 4, 3"), awk = NULL)
    awkHelper(file, "gear contains 3", awk = NULL)
    # note that to be filtered the file has to exist on your local system, this example only shows
    # the output of awkHelper, which would then be executed by read.pcv on a unix system
    awkHelper(file, list("gear in 4, 3"), awk = "existing_command")
  },
  error = function(e) {
    message(e)
  }
)

Description

The Bayesian Analysis Reporting Guidelines were put forward by Kruschke (https://www.nature.com/articles/s41562-021-01177-7) to aide in reproducibility and documentation of Bayesian statistical analyses that are sometimes unfamiliar to reviewers or scientists. The purpose of this function is to summarize goodness of fit metrics from one or more Bayesian models made by growthSS and fitGrowth. See details for explanations of those metrics and the output.

Usage

barg(fit, ss = NULL)

Arguments

Details

• General: This includes chain number, length, and total divergent transitions per model. Divergent transitions are a marker that the MCMC had something go wrong. Conceptually it may be helpful to think about rolling a marble over a 3D curve then having the marble suddenly jolt in an unexpected direction, something happened that suggests a problem/misunderstood surface. In practice you want extremely few (ideally no) divergences. If you do have divergences then consider specifying more control parameters (see brms::brm or examples for fitGrowth). If the problem persists then the model may need to be simplified. For more information on MCMC and divergence see the stan manual (https://mc-stan.org/docs/2_19/reference-manual/divergent-transitions).

• ESS: ESS stands for Effective Sample Size and is a goodness of fit metric that approximates the number of independent replicates that would equate to the same amount of information as the (autocorrelated) MCMC iterations. ESS of 1000+ is often considered as a pretty stable value, but more is better. Still, 100 per chain may be plenty depending on your applications and the inference you wish to do. One of the benefits to using lots of chains and/or longer chains is that you will get more complete information and that benefit will be shown by a larger ESS. This is separated into "bulk" and "tail" to represent the middle and tails of the posterior distribution, since those can sometimes have very different sampling behavior. A summary and the total values are returned, with the summary being useful if several models are included in a list for fit argument

• Rhat: Rhat is a measure of "chain mixture". It compares the between vs within chain values to assess how well the chains mixed. If chains did not mix well then Rhat will be greater than 1, with 1.05 being a broadly agreed upon cutoff to signify a problem. Running longer chains should result in lower Rhat values. The default in brms is to run 4 chains, partially to ensure that there is a good chance to check that the chains mixed well via Rhat. A summary and the total values are returned, with the summary being useful if several models are included in a list for fit argument

• NEFF: NEFF is the NEFF ratio (Effective Sample Size over Total MCMC Sample Size). Values greater than 0.5 are generally considered good, but there is a consensus that lower can be fine down to about 0.1. A summary and the total values are returned, with the summary being useful if several models are included in a list for fit argument

• priorPredictive: A plot of data simulated from the prior using plotPrior. This should generate data that is biologically plausible for your situation, but it will probably be much more variable than your data. That is the effect of the mildly informative thick tailed lognormal priors. If you specified non-default style priors then this currently will not work.

• posteriorPredictive: A plot of each model's posterior predictive interval over time. This is the same as plots returned from growthPlot and shows 1-99 coming to a mean yellow trend line. These should encompass the overwhelming majority of your data and ideally match the variance pattern that you see in your data. If parts of the predicted interval are biologically impossible (area below 0, percentage about 100 model should be reconsidered.

Value

A named list containing Rhat, ESS, NEFF, and Prior/Posterior Predictive plots. See details for interpretation.

See Also

plotPrior for visual prior predictive checks.

Examples

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
)
barg(fit_test, ss)
fit_2 <- fit_test
fit_list <- list(fit_test, fit_2)
x <- barg(fit_list, list(ss, ss))

Description

Models fit using growthSS inputs by fitGrowth (and similar models made through other means) can be visualized easily using this function. This will generally be called by growthPlot.

Usage

brmPlot(
  fit,
  form,
  df = NULL,
  groups = NULL,
  timeRange = NULL,
  facetGroups = TRUE,
  hierarchy_value = NULL,
  vir_option = "plasma"
)

Arguments

Value

Returns a ggplot showing a brms model's credible intervals and optionally the individual growth lines.

Examples

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)
growthPlot(fit = fit, form = y ~ time | group, groups = "a", df = ss$df)

Description

Models fit using growthSS inputs by fitGrowth (and similar models made through other means) can be visualized easily using this function. This will generally be called by growthPlot.

Usage

brmSurvPlot(
  fit,
  form,
  df = NULL,
  groups = NULL,
  timeRange = NULL,
  facetGroups = TRUE
)

Arguments

Value

Returns a ggplot showing a brms model's credible intervals and optionally the individual growth lines.

Examples

set.seed(123)
df <- growthSim("exponential",
  n = 20, t = 50,
  params = list("A" = c(1, 1), "B" = c(0.15, 0.1))
)
ss1 <- 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)

Description

Function to visualize hypotheses tested on brms models similar to those made using growthSS outputs.

Usage

brmViolin(fit, ss, hypothesis)

Arguments

Value

Returns a ggplot showing a brms model's posterior distributions as violins and filled by posterior probability of some hypothesis.

Examples

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)
brmViolin(fit, ss, ".../A_groupd > 1.05") # all groups used
brmViolin(fit, ss, "A_groupa/A_groupd > 1.05") # only these two groups

Description

Remove outliers from bellwether data using cook's distance

Usage

bw.outliers(
  df = NULL,
  phenotype,
  naTo0 = FALSE,
  group = c(),
  cutoff = 3,
  outlierMethod = "cooks",
  plotgroup = c("barcode", "rotation"),
  plot = TRUE,
  x = NULL,
  traitCol = "trait",
  valueCol = "value",
  labelCol = "label",
  idCol = NULL,
  ncp = NULL,
  separate = NULL
)

Arguments

Value

The input dataframe with outliers removed and optionally a plot (if a plot is returned then output is a list).

Examples

sv <- growthSim("logistic",
  n = 5, t = 20,
  params = list("A" = c(200, 160), "B" = c(13, 11), "C" = c(3, 3.5))
)
sv[130, ]$y <- 500
sv_res <- bw.outliers(
  df = sv, phenotype = "y", naTo0 = FALSE, cutoff = 15,
  group = c("time", "group"), outlierMethod = "cooks",
  plotgroup = "id", plot = TRUE
)
sv_res$plot

tryCatch(
  { # in case offline
    library(data.table)
    mvw <- read.pcv(paste0(
      "https://media.githubusercontent.com/media/joshqsumner/",
      "pcvrTestData/main/pcv4-multi-value-traits.csv"
    ), mode = "wide", reader = "fread")
    mvw$genotype <- substr(mvw$barcode, 3, 5)
    mvw$genotype <- ifelse(mvw$genotype == "002", "B73",
      ifelse(mvw$genotype == "003", "W605S",
        ifelse(mvw$genotype == "004", "MM", "Mo17")
      )
    )
    mvw$fertilizer <- substr(mvw$barcode, 8, 8)
    mvw$fertilizer <- ifelse(mvw$fertilizer == "A", "100",
      ifelse(mvw$fertilizer == "B", "50", "0")
    )
    mvw <- bw.time(mvw, timeCol = "timestamp", group = "barcode", plot = FALSE)

    phenotypes <- which(grepl("hue_freq", colnames(mvw)))

    mvw2 <- bw.outliers(
      df = mvw, phenotype = phenotypes, naTo0 = FALSE, outlierMethod = "cooks",
      group = c("DAS", "genotype", "fertilizer"), cutoff = 3, plotgroup = c("barcode", "rotation")
    )


    mvl <- read.pcv(paste0(
      "https://media.githubusercontent.com/media/joshqsumner/",
      "pcvrTestData/main/pcv4-multi-value-traits.csv"
    ), mode = "long")
    mvl$genotype <- substr(mvl$barcode, 3, 5)
    mvl$genotype <- ifelse(mvl$genotype == "002", "B73",
      ifelse(mvl$genotype == "003", "W605S",
        ifelse(mvl$genotype == "004", "MM", "Mo17")
      )
    )
    mvl$fertilizer <- substr(mvl$barcode, 8, 8)
    mvl$fertilizer <- ifelse(mvl$fertilizer == "A", "100",
      ifelse(mvl$fertilizer == "B", "50", "0")
    )
    mvl <- bw.time(mvl, timeCol = "timestamp", group = "barcode", plot = FALSE)

    mvl2 <- bw.outliers(
      df = mvl, phenotype = "hue_frequencies", naTo0 = FALSE, outlierMethod = "cooks",
      group = c("DAS", "genotype", "fertilizer"), cutoff = 3, plotgroup = c("barcode", "rotation")
    )
  },
  error = function(e) {
    message(e)
  }
)

Description

Time conversion and plotting for bellwether data

Usage

bw.time(
  df = NULL,
  mode = c("DAS", "DAP", "DAE"),
  plantingDelay = NULL,
  phenotype = NULL,
  cutoff = 1,
  timeCol = "timestamp",
  group = "Barcodes",
  plot = TRUE,
  format = "%Y-%m-%d %H:%M:%S",
  traitCol = "trait",
  valueCol = "value",
  index = NULL
)

Arguments

Value

The input dataframe with new integer columns for different ways of describing time in the experiment. If plot is TRUE then a ggplot is also returned as part of a list.

Examples

f <- "https://raw.githubusercontent.com/joshqsumner/pcvrTestData/main/pcv4-single-value-traits.csv"
tryCatch(
  {
    sv <- read.pcv(
      f,
      mode = "wide", reader = "fread"
    )
    sv$genotype = substr(sv$barcode, 3, 5)
    sv$genotype = ifelse(sv$genotype == "002", "B73",
      ifelse(sv$genotype == "003", "W605S",
        ifelse(sv$genotype == "004", "MM", "Mo17")
      )
    )
    sv$fertilizer = substr(sv$barcode, 8, 8)
    sv$fertilizer = ifelse(sv$fertilizer == "A", "100",
      ifelse(sv$fertilizer == "B", "50", "0")
    )
    sv <- bw.time(sv,
      plantingDelay = 0, phenotype = "area_pixels", cutoff = 10,
      timeCol = "timestamp", group = c("barcode", "rotation"), plot = FALSE
    )

    svl <- read.pcv(
      f,
      mode = "long", reader = "fread"
    )
    svl$genotype = substr(svl$barcode, 3, 5)
    svl$genotype = ifelse(svl$genotype == "002", "B73",
      ifelse(svl$genotype == "003", "W605S",
        ifelse(svl$genotype == "004", "MM", "Mo17")
      )
    )
    svl$fertilizer = substr(svl$barcode, 8, 8)
    svl$fertilizer = ifelse(svl$fertilizer == "A", "100",
      ifelse(svl$fertilizer == "B", "50", "0")
    )
    svl <- bw.time(svl,
      plantingDelay = 0, phenotype = "area_pixels", cutoff = 10, timeCol = "timestamp",
      group = c("barcode", "rotation"), plot = FALSE
    )
  },
  error = function(e) {
    message(e)
  }
)

Description

Read in lemnatech watering data from metadata.json files

Usage

bw.water(file = NULL, envKey = "environment")

Arguments

Value

A data frame containing the bellwether watering data

Examples

tryCatch(
  {
    w <- bw.water("https://raw.githubusercontent.com/joshqsumner/pcvrTestData/main/metadata.json")
  },
  error = function(e) {
    message(e)
  }
)

Description

Helper function to check groups in data.

Usage

checkGroups(df, group)

Arguments

Value

If there are duplicates in the grouping then this will return a message with code to start checking the duplicates in your data.

Examples

df <- growthSim("linear",
  n = 10, t = 10,
  params = list("A" = c(2, 1.5))
)
checkGroups(df, c("time", "id", "group"))
df$time[12] <- 3
checkGroups(df, c("time", "id", "group"))

Description

Helper function for binding draws from several brms models to make a data.frame for use with brms::hypothesis(). This will also check that the draws are comparable using basic model metrics.

Usage

combineDraws(..., message = TRUE)

Arguments

Value

Returns a dataframe of posterior draws.

Examples

# note that this example will fit several bayesian models and may run for several minutes.

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)

Description

Function to perform bayesian tests and ROPE comparisons using single or multi value traits with several distributions.

Usage

conjugate(
  s1 = NULL,
  s2 = NULL,
  method = c("t", "gaussian", "beta", "binomial", "lognormal", "lognormal2", "poisson",
    "negbin", "vonmises", "vonmises2", "uniform", "pareto", "gamma", "bernoulli",
    "exponential", "bivariate_uniform", "bivariate_gaussian", "bivariate_lognormal"),
  priors = NULL,
  plot = FALSE,
  rope_range = NULL,
  rope_ci = 0.89,
  cred.int.level = 0.89,
  hypothesis = "equal",
  support = NULL
)

Arguments

Details

Prior distributions default to be weakly informative and in some cases you may wish to change them.

• "t" and "gaussian": priors = list( mu=c(0,0),n=c(1,1),s2=c(20,20) ) , where mu is the mean, n is the number of prior observations, and s2 is variance

• "beta", "bernoulli", and "binomial": priors = list( a=c(0.5, 0.5), b=c(0.5, 0.5) ), where a and b are shape parameters of the beta distribution. Note that for the binomial distribution this is used as the prior for success probability P, which is assumed to be beta distributed as in a beta-binomial distribution.

• "lognormal": priors = list(mu = 0, sd = 5) , where mu and sd describe the normal distribution of the mean parameter for lognormal data. Note that these values are on the log scale.

• "lognormal2": priors = list(a = 1, b = 1) , where a and b are the shape and scale parameters of the gamma distribution of lognormal data's precision parameter (using the alternative mu, precision paramterization).

• "gamma": priors = list(shape = 0.5, scale = 0.5, known_shape = 1), where shape and scale are the respective parameters of the gamma distributed rate (inverse of scale) parameter of gamma distributed data.

• "poisson" and "exponential": priors = list(a=c(0.5,0.5),b=c(0.5,0.5)), where a and b are shape parameters of the gamma distribution.

• "negbin": priors = list(r=c(10,10), a=c(0.5,0.5),b=c(0.5,0.5)), where r is the r parameter of the negative binomial distribution (representing the number of successes required) and where a and b are shape parameters of the beta distribution. Note that the r value is not updated. The conjugate beta prior is only valid when r is fixed and known, which is a limitation for this method.

• "uniform": list(scale = 0.5, location = 0.5), where scale is the scale parameter of the pareto distributed upper boundary and location is the location parameter of the pareto distributed upper boundary. Note that different sources will use different terminology for these parameters. These names were chosen for consistency with the extraDistr implementation of the pareto distribution. On Wikipedia the parameters are called shape and scale, corresponding to extraDistr's scale and location respecitvely, which can be confusing. Note that the lower boundary of the uniform is assumed to be 0.

• "pareto": list(a = 1, b = 1, known_location = min(data)), where a and b are the shape and scale parameters of the gamma distribution of the pareto distribution's scale parameter. In this case location is assumed to be constant and known, which is less of a limitation than knowing r for the negative binomial method since location will generally be right around/just under the minimum of the sample data. Note that the pareto method is only implemented currently for single value traits since one of the statistics needed to update the gamma distribution here is the product of the data and we do not currently have a method to calculate a similar sufficient statistic from multi value traits.

• "vonmises": list(mu = 0, kappa = 0.5, boundary = c(-pi, pi), known_kappa = 1, n = 1), where mu is the direction of the circular distribution (the mean), kappa is the precision of the mean, boundary is a vector including the two values that are the where the circular data "wraps" around the circle, known_kappa is the fixed value of precision for the total distribution, and n is the number of prior observations. This Von Mises option updates the conjugate prior for the mean direction, which is itself Von-Mises distributed. This in some ways is analogous to the T method, but assuming a fixed variance when the mean is updated. Note that due to how the rescaling works larger circular boundaries can be slow to plot.

• "vonmises2": priors = list(mu = 0, kappa = 0.5, boundary = c(-pi, pi), n = 1), where mu and kappa are mean direction and precision of the von mises distribution, boundary is a vector including the two values that are the where the circular data "wraps" around the circle, and n is the number of prior observations. This Von-Mises implementation does not assume constant variance and instead uses MLE to estimate kappa from the data and updates the kappa prior as a weighted average of the data and the prior. The mu parameter is then updated per Von-Mises conjugacy.

• "bivariate_uniform": list(location_l = 1, location_u = 2, scale = 1), where scale is the shared scale parameter of the pareto distributed upper and lower boundaries and location l and u are the location parameters for the Lower (l) and Upper (u) boundaries of the uniform distribution. Note this uses the same terminology for the pareto distribution's parameters as the "uniform" method.

• "bivariate_gaussian" and "bivariate_lognormal": list(mu = 0, sd = 10, a = 1, b = 1), where mu and sd are the mean and standard deviation of the Normal distribution of the data's mean and a and b are the shape and scale of the gamma distribution on precision. Note that internally this uses the Mu and Precision parameterization of the normal distribution and those are the parameters shown in the plot and tested, but priors use Mu and SD for the normal distribution of the mean.

See examples for plots of these prior distributions.

Value

A list with named elements:

• summary: A data frame containing HDI/HDE values for each sample and the ROPE as well as posterior probability of the hypothesis and ROPE test (if specified). The HDE is the "Highest Density Estimate" of the posterior, that is the tallest part of the probability density function. The HDI is the Highest Density Interval, which is an interval that contains X% of the posterior distribution, so cred.int.level = 0.8 corresponds to an HDI that includes 80 percent of the posterior probability.

• posterior: A list of updated parameters in the same format as the prior for the given method. If desired this does allow for Bayesian updating.

• plot_df: A data frame of probabilities along the support for each sample. This is used for making the ggplot.

• rope_df: A data frame of draws from the ROPE posterior.

• plot: A ggplot showing the distribution of samples and optionally the distribution of differences/ROPE

Examples

mv_ln <- mvSim(
  dists = list(
    rlnorm = list(meanlog = log(130), sdlog = log(1.2)),
    rlnorm = list(meanlog = log(100), sdlog = log(1.3))
  ),
  n_samples = 30
)

# lognormal mv
ln_mv_ex <- conjugate(
  s1 = mv_ln[1:30, -1], s2 = mv_ln[31:60, -1], method = "lognormal",
  priors = list(mu = 5, sd = 2),
  plot = FALSE, rope_range = c(-40, 40), rope_ci = 0.89,
  cred.int.level = 0.89, hypothesis = "equal", support = NULL
)

# lognormal sv
ln_sv_ex <- conjugate(
  s1 = rlnorm(100, log(130), log(1.3)), s2 = rlnorm(100, log(100), log(1.6)),
  method = "lognormal",
  priors = list(mu = 5, sd = 2),
  plot = FALSE, rope_range = NULL, rope_ci = 0.89,
  cred.int.level = 0.89, hypothesis = "equal", support = NULL
)

# Z test mv example

mv_gauss <- mvSim(
  dists = list(
    rnorm = list(mean = 50, sd = 10),
    rnorm = list(mean = 60, sd = 12)
  ),
  n_samples = 30
)

gauss_mv_ex <- conjugate(
  s1 = mv_gauss[1:30, -1], s2 = mv_gauss[31:60, -1], method = "gaussian",
  priors = list(mu = 30, n = 1, s2 = 100),
  plot = FALSE, rope_range = c(-25, 25), rope_ci = 0.89,
  cred.int.level = 0.89, hypothesis = "equal", support = NULL
)

# T test sv example

gaussianMeans_sv_ex <- conjugate(
  s1 = rnorm(10, 50, 10), s2 = rnorm(10, 60, 12), method = "t",
  priors = list(mu = 30, n = 1, s2 = 100),
  plot = FALSE, rope_range = c(-5, 8), rope_ci = 0.89,
  cred.int.level = 0.89, hypothesis = "equal", support = NULL
)

# beta mv example

set.seed(123)
mv_beta <- mvSim(
  dists = list(
    rbeta = list(shape1 = 5, shape2 = 8),
    rbeta = list(shape1 = 10, shape2 = 10)
  ),
  n_samples = c(30, 20)
)

beta_mv_ex <- conjugate(
  s1 = mv_beta[1:30, -1], s2 = mv_beta[31:50, -1], method = "beta",
  priors = list(a = 0.5, b = 0.5),
  plot = FALSE, rope_range = c(-0.1, 0.1), rope_ci = 0.89,
  cred.int.level = 0.89, hypothesis = "equal"
)

# beta sv example

beta_sv_ex <- conjugate(
  s1 = rbeta(20, 5, 5), s2 = rbeta(20, 8, 5), method = "beta",
  priors = list(a = 0.5, b = 0.5),
  plot = FALSE, rope_range = c(-0.1, 0.1), rope_ci = 0.89,
  cred.int.level = 0.89, hypothesis = "equal"
)

# binomial sv example
# note that specifying trials = 20 would also work
# and the number of trials will be recycled to the length of successes

binomial_sv_ex <- conjugate(
  s1 = list(successes = c(15, 14, 16, 11), trials = c(20, 20, 20, 20)),
  s2 = list(successes = c(7, 8, 10, 5), trials = c(20, 20, 20, 20)), method = "binomial",
  priors = list(a = 0.5, b = 0.5),
  plot = FALSE, rope_range = c(-0.1, 0.1), rope_ci = 0.89,
  cred.int.level = 0.89, hypothesis = "equal"
)

# poisson sv example

poisson_sv_ex <- conjugate(
  s1 = rpois(20, 10), s2 = rpois(20, 8), method = "poisson",
  priors = list(a = 0.5, b = 0.5),
  plot = FALSE, rope_range = c(-1, 1), rope_ci = 0.89,
  cred.int.level = 0.89, hypothesis = "equal"
)

# negative binomial sv example
# knowing r (required number of successes) is an important caveat for this method.
# in the current implementation we suggest using the poisson method for data such as leaf counts

negbin_sv_ex <- conjugate(
  s1 = rnbinom(20, 10, 0.5), s2 = rnbinom(20, 10, 0.25), method = "negbin",
  priors = list(r = 10, a = 0.5, b = 0.5),
  plot = FALSE, rope_range = c(-1, 1), rope_ci = 0.89,
  cred.int.level = 0.89, hypothesis = "equal"
)

# von mises mv example

mv_gauss <- mvSim(
  dists = list(
    rnorm = list(mean = 50, sd = 10),
    rnorm = list(mean = 60, sd = 12)
  ),
  n_samples = c(30, 40)
)
vm1_ex <- conjugate(
  s1 = mv_gauss[1:30, -1],
  s2 = mv_gauss[31:70, -1],
  method = "vonmises",
  priors = list(mu = 45, kappa = 1, boundary = c(0, 180), known_kappa = 1, n = 1),
  plot = FALSE, rope_range = c(-1, 1), rope_ci = 0.89,
  cred.int.level = 0.89, hypothesis = "equal"
)

# von mises 2 sv example
vm2_ex <- conjugate(
  s1 = brms::rvon_mises(10, 2, 2),
  s2 = brms::rvon_mises(15, 3, 3),
  method = "vonmises2",
  priors = list(mu = 0, kappa = 0.5, boundary = c(-pi, pi), n = 1),
  cred.int.level = 0.95,
  plot = FALSE
)

Description

Often in bellwether experiments we are curious about the effect of some treatment vs control. For certain routes in analysing the data this requires considering phenotypes as relative differences compared to a control.

Usage

cumulativePheno(
  df,
  phenotypes = NULL,
  group = "barcode",
  timeCol = "DAS",
  traitCol = "trait",
  valueCol = "value"
)

Arguments

Value

A dataframe with cumulative sum columns added for specified phenotypes

Examples

f <- "https://raw.githubusercontent.com/joshqsumner/pcvrTestData/main/pcv4-single-value-traits.csv"
tryCatch(
  {
    sv <- read.pcv(
      f,
      reader = "fread"
    )
    sv$genotype <- substr(sv$barcode, 3, 5)
    sv$genotype <- ifelse(sv$genotype == "002", "B73",
      ifelse(sv$genotype == "003", "W605S",
        ifelse(sv$genotype == "004", "MM", "Mo17")
      )
    )
    sv$fertilizer <- substr(sv$barcode, 8, 8)
    sv$fertilizer <- ifelse(sv$fertilizer == "A", "100",
      ifelse(sv$fertilizer == "B", "50", "0")
    )

    sv <- bw.time(sv,
      plantingDelay = 0, phenotype = "area_pixels", cutoff = 10,
      timeCol = "timestamp", group = c("barcode", "rotation"), plot = TRUE
    )$data
    sv <- bw.outliers(sv,
      phenotype = "area_pixels", group = c("DAS", "genotype", "fertilizer"),
      plotgroup = c("barcode", "rotation")
    )$data
    phenotypes <- colnames(sv)[19:35]
    phenoForm <- paste0("cbind(", paste0(phenotypes, collapse = ", "), ")")
    groupForm <- "DAS+DAP+barcode+genotype+fertilizer"
    form <- as.formula(paste0(phenoForm, "~", groupForm))
    sv <- aggregate(form, data = sv, mean, na.rm = TRUE)
    pixels_per_cmsq <- 42.5^2 # pixel per cm^2
    sv$area_cm2 <- sv$area_pixels / pixels_per_cmsq
    sv$height_cm <- sv$height_pixels / 42.5
    df <- sv
    phenotypes <- c("area_cm2", "height_cm")
    group <- c("barcode")
    timeCol <- "DAS"
    df <- cumulativePheno(df, phenotypes, group, timeCol)


    sv_l <- read.pcv(
      f,
      mode = "long", reader = "fread"
    )
    sv_l$genotype <- substr(sv_l$barcode, 3, 5)
    sv_l$genotype <- ifelse(sv_l$genotype == "002", "B73",
      ifelse(sv_l$genotype == "003", "W605S",
        ifelse(sv_l$genotype == "004", "MM", "Mo17")
      )
    )
    sv_l$fertilizer <- substr(sv_l$barcode, 8, 8)
    sv_l$fertilizer <- ifelse(sv_l$fertilizer == "A", "100",
      ifelse(sv_l$fertilizer == "B", "50", "0")
    )
    sv_l <- bw.time(sv_l,
      plantingDelay = 0, phenotype = "area_pixels", cutoff = 10,
      timeCol = "timestamp", group = c("barcode", "rotation")
    )$data
    sv_l <- cumulativePheno(sv_l,
      phenotypes = c("area_pixels", "height_pixels"),
      group = c("barcode", "rotation"), timeCol = "DAS"
    )
  },
  error = function(e) {
    message(e)
  }
)

Description

Function for plotting iterations of posterior distributions

Usage

distributionPlot(
  fits,
  form,
  df,
  priors = NULL,
  params = NULL,
  maxTime = NULL,
  patch = TRUE
)

Arguments

Value

A ggplot or a list of ggplots (depending on patch).

Examples

f <- "https://raw.githubusercontent.com/joshqsumner/pcvrTestData/main/brmsFits.rdata"
tryCatch(
  {
    print(load(url(f)))
    library(brms)
    library(ggplot2)
    library(patchwork)
    fits <- list(fit_3, fit_15)
    form <- y~time | id / group
    priors <- list(
      "phi1" = rlnorm(2000, log(130), 0.25),
      "phi2" = rlnorm(2000, log(12), 0.25),
      "phi3" = rlnorm(2000, log(3), 0.25)
    )
    params <- c("A", "B", "C")
    d <- simdf
    maxTime <- NULL
    patch <- TRUE
    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 = params, d = d, priors = priors, patch = FALSE
    )
    distributionPlot(
      fits = from3to25, form = y ~ time | id / group,
      params = params, d = d, patch = FALSE
    )
  },
  error = function(e) {
    message(e)
  }
)

## End(Not run)

Ease of use wrapper function for fitting various growth models specified by growthSS

Description

Ease of use wrapper function for fitting various growth models specified by growthSS

Usage

fitGrowth(ss, ...)

Arguments

Value

A fit model from the selected type.

See Also

growthPlot for model visualization, testGrowth for hypothesis testing, barg for Bayesian model reporting metrics.

Examples

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 | group,
  df = simdf, type = "nls"
)
fitGrowth(ss)
ss <- growthSS(
  model = "gam", form = y ~ time | group,
  df = simdf, type = "nls"
)
fitGrowth(ss)

Description

Helper function generally called from fitGrowth.

Usage

fitGrowthbrms(
  ss,
  iter = 2000,
  cores = getOption("mc.cores", 1),
  chains = 4,
  prior = NULL,
  backend = "cmdstanr",
  silent = 0,
  ...
)

fitGrowthbrmsgam(
  ss,
  iter = 2000,
  cores = getOption("mc.cores", 1),
  chains = 4,
  prior = NULL,
  backend = "cmdstanr",
  silent = 0,
  ...
)

Arguments

Value

A brmsfit object, see ?`brmsfit-class` for details.

Description

Helper function generally called from fitGrowth.

Usage

fitGrowthflexsurv(ss, ...)

Arguments

Value

A survreg object.

Description

Helper function generally called from fitGrowth.

Usage

fitGrowthlm(ss, ...)

Arguments

Value

An lm object, see ?lm for details.

Description

Helper function generally called from fitGrowth.

Usage

fitGrowthmgcvgam(ss, ...)

Arguments

Value

An gam object, see ?gam for details.

Description

Helper function generally called from fitGrowth.

Usage

fitGrowthnlme(ss, ...)

Arguments

Value

An nlme object, see ?nlme for details.

Description

Helper function generally called from fitGrowth.

Usage

fitGrowthnlmegam(ss, ...)

Arguments

Value

An lme object, see ?lme for details.

Description

Helper function generally called from fitGrowth.

Usage

fitGrowthnlrq(ss, cores = getOption("mc.cores", 1), ...)

Arguments

Value

An nlrqModel object if tau is length of 1 or a list of such models for multiple taus, see ?quantreg::nlrq or ?nls::nlsModel for details.

Description

Helper function generally called from fitGrowth.

Usage

fitGrowthnlrqgam(ss, cores = getOption("mc.cores", 1), ...)

Arguments

Value

An rq object, see ?rq for details.

Description

Helper function generally called from fitGrowth.

Usage

fitGrowthnls(ss, ...)

Arguments

Value

An nls object, see ?nls for details.

Description

Helper function generally called from fitGrowth.

Usage

fitGrowthnlsgam(ss, ...)

Arguments

Value

An lm object, see ?lm for details.

Description

Helper function generally called from fitGrowth.

Usage

fitGrowthrq(ss, cores = getOption("mc.cores", 1), ...)

Arguments

Value

An rq object, see ?rq for details.

Description

Helper function generally called from fitGrowth.

Usage

fitGrowthsurvreg(ss, ...)

Arguments

Value

A survreg object.

Description

Models fit using growthSS inputs by fitGrowth (and similar models made through other means) can be visualized easily using this function. This will generally be called by growthPlot.

Usage

flexsurvregPlot(
  fit,
  form,
  groups = NULL,
  df = NULL,
  timeRange = NULL,
  facetGroups = TRUE,
  groupFill = FALSE,
  virMaps = c("plasma")
)

Arguments

Value

Returns a ggplot showing an survival model's survival function.

Examples

df <- growthSim("logistic",
  n = 20, t = 25,
  params = list("A" = c(200, 160), "B" = c(13, 11), "C" = c(3, 3.5))
)
ss <- growthSS(
  model = "survival weibull", form = y > 100 ~ time | id / group,
  df = df, type = "flexsurv"
)
fit <- fitGrowth(ss)
flexsurvregPlot(fit, form = ss$pcvrForm, df = ss$df, groups = "a")
flexsurvregPlot(fit,
  form = ss$pcvrForm, df = ss$df,
  facetGroups = FALSE, groupFill = TRUE
)

Description

Variance partitioning for phenotypes (over time) using fully random effects models

Usage

frem(
  df,
  des,
  phenotypes,
  timeCol = NULL,
  cor = TRUE,
  returnData = FALSE,
  combine = TRUE,
  markSingular = FALSE,
  time = NULL,
  time_format = "%Y-%m-%d",
  ...
)

Arguments

Value

Returns either a plot (if returnData=FALSE) or a list with a plot and data/a list of dataframes (depending on returnData and cor).

Examples

library(data.table)
set.seed(456)
df <- data.frame(
  genotype = rep(c("g1", "g2"), each = 10),
  treatment = rep(c("C", "T"), times = 10),
  time = rep(c(1:5), times = 2),
  date_time = rep(paste0("2024-08-", 21:25), times = 2),
  pheno1 = rnorm(20, 10, 1),
  pheno2 = sort(rnorm(20, 5, 1)),
  pheno3 = sort(runif(20))
)
out <- frem(df, des = "genotype", phenotypes = c("pheno1", "pheno2", "pheno3"), returnData = TRUE)
lapply(out, class)
frem(df,
  des = c("genotype", "treatment"), phenotypes = c("pheno1", "pheno2", "pheno3"),
  cor = FALSE
)
frem(df,
  des = "genotype", phenotypes = c("pheno1", "pheno2", "pheno3"),
  combine = FALSE, timeCol = "time", time = "all"
)
frem(df,
  des = "genotype", phenotypes = c("pheno1", "pheno2", "pheno3"),
  combine = TRUE, timeCol = "time", time = 1
)
frem(df,
  des = "genotype", phenotypes = c("pheno1", "pheno2", "pheno3"),
  cor = FALSE, timeCol = "time", time = 3:5, markSingular = TRUE
)
df[df$time == 3, "genotype"] <- "g1"
frem(df,
  des = "genotype", phenotypes = c("pheno1", "pheno2", "pheno3"),
  cor = FALSE, timeCol = "date_time", time = "all", markSingular = TRUE
)

Description

Note that using GAMs will be less useful than fitting parameterized models as supported by growthSS and fitGrowth for common applications in plant phenotyping.

Usage

gam_diff(
  model,
  newdata = NULL,
  g1,
  g2,
  byVar = NULL,
  smoothVar = NULL,
  cis = seq(0.05, 0.95, 0.05),
  unconditional = TRUE,
  plot = TRUE
)

Arguments

Value

A dataframe or a list containing a ggplot and a dataframe

Examples

ex <- pcvr::growthSim("logistic",
  n = 20, t = 25,
  params = list(
    "A" = c(200, 160),
    "B" = c(13, 11),
    "C" = c(3, 3.5)
  )
)

m <- mgcv::gam(y ~ group + s(time, by = factor(group)), data = ex)

support <- expand.grid(
  time = seq(min(ex$time), max(ex$time), length = 400),
  group = factor(unique(ex$group))
)

out <- gam_diff(
  model = m, newdata = support, g1 = "a", g2 = "b",
  byVar = "group", smoothVar = "time", plot = TRUE
)
dim(out$data)
out$plot
out2 <- gam_diff(
  model = m, g1 = "a", g2 = "b", byVar = NULL, smoothVar = NULL, plot = TRUE
)

Function to visualize models made by fitGrowth.

Description

Models fit using growthSS inputs by fitGrowth (and similar models made through other means) can be visualized easily using this function.

Usage

growthPlot(
  fit,
  form,
  groups = NULL,
  df = NULL,
  timeRange = NULL,
  facetGroups = TRUE,
  groupFill = !facetGroups,
  hierarchy_value = NULL
)

Arguments

Value

Returns a ggplot showing a brms model's credible intervals and optionally the individual growth lines.

See Also

growthSS and fitGrowth for making compatible models, testGrowth for hypothesis testing on compatible models.

Examples

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, type = "nls"
)
fit <- fitGrowth(ss)
growthPlot(fit, form = ss$pcvrForm, df = ss$df)

Description

growthSim can be used to help pick reasonable parameters for common growth models to use in prior distributions or to simulate data for example models/plots.

Usage

growthSim(
  model = c("logistic", "gompertz", "double logistic", "double gompertz",
    "monomolecular", "exponential", "linear", "power law", "frechet", "weibull",
    "gumbel", "logarithmic", "bragg", "lorentz", "beta"),
  n = 20,
  t = 25,
  params = list(),
  D = 0
)

Arguments

Details

The params argument requires some understanding of how each growth model is parameterized. Examples of each are below should help, as will the examples.

• Logistic: 'A / (1 + exp( (B-x)/C) )' Where A is the asymptote, B is the inflection point, C is the growth rate.

• Gompertz: 'A * exp(-B * exp(-C*x))' Where A is the asymptote, B is the inflection point, C is the growth rate.

• Weibull: 'A * (1-exp(-(x/C)^B))' Where A is the asymptote, B is the weibull shape parameter, C is the weibull scale parameter.

• Frechet: 'A * exp(-((x-0)/C)^(-B))' Where A is the asymptote, B is the frechet shape parameter, C is the frechet scale parameter. Note that the location parameter (conventionally m) is 0 in these models for simplicity but is still included in the formula.

• Gumbel: 'A * exp(-exp(-(x-B)/C))' Where A is the asymptote, B is the inflection point (location), C is the growth rate (scale).

• Double Logistic: 'A / (1+exp((B-x)/C)) + ((A2-A) /(1+exp((B2-x)/C2)))' Where A is the asymptote, B is the inflection point, C is the growth rate, A2 is the second asymptote, B2 is the second inflection point, and C2 is the second growth rate.

• Double Gompertz: 'A * exp(-B * exp(-C*x)) + ((A2-A) * exp(-B2 * exp(-C2*(x-B))))' Where A is the asymptote, B is the inflection point, C is the growth rate, A2 is the second asymptote, B2 is the second inflection point, and C2 is the second growth rate.

• Monomolecular: 'A-A * exp(-B * x)' Where A is the asymptote and B is the growth rate.

• Exponential: 'A * exp(B * x)' Where A is the scale parameter and B is the growth rate.

• Linear: 'A * x' Where A is the growth rate.

• Logarithmic: 'A * log(x)' Where A is the growth rate.

• Power Law: 'A * x ^ (B)' Where A is the scale parameter and B is the growth rate.

• Bragg: 'A * exp(-B * (x - C) ^ 2)' This models minima and maxima as a dose-response curve where A is the max response, B is the "precision" or slope at inflection, and C is the x position of the max response.

• Lorentz: 'A / (1 + B * (x - C) ^ 2)' This models minima and maxima as a dose-response curve where A is the max response, B is the "precision" or slope at inflection, and C is the x position of the max response. Generally Bragg is preferred to Lorentz for dose-response curves.

• Beta: 'A * (((x - D) / (C - D)) * ((E - x) / (E - C)) ^ ((E - C) / (C - D))) ^ B' This models minima and maxima as a dose-response curve where A is the Maximum Value, B is a shape/concavity exponent similar to the sum of alpha and beta in a Beta distribution, C is the position of maximum value, D is the minimum position where distribution > 0, E is the maximum position where distribution > 0. This is a difficult model to fit but can model non-symmetric dose-response relationships which may sometimes be worth the extra effort.

Note that for these distributions parameters generally do not exist in a vacuum. Changing one will make the others look different in the resulting data. The examples are a good place to start if you are unsure what parameters to use.

Value

Returns a dataframe of example growth data following the input parameters.

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))
)
ggplot(simdf, aes(time, y, group = interaction(group, id))) +
  geom_line(aes(color = group)) +
  labs(title = "Logistic")

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

simdf <- growthSim("weibull",
  n = 20, t = 25,
  params = list("A" = c(100, 100), "B" = c(1, 0.75), "C" = c(2, 3))
)
ggplot(simdf, aes(time, y, group = interaction(group, id))) +
  geom_line(aes(color = group)) +
  labs(title = "weibull")

simdf <- growthSim("frechet",
  n = 20, t = 25,
  params = list("A" = c(100, 110), "B" = c(2, 1.5), "C" = c(5, 2))
)
ggplot(simdf, aes(time, y, group = interaction(group, id))) +
  geom_line(aes(color = group)) +
  labs(title = "frechet")

simdf <- growthSim("gumbel",
  n = 20, t = 25,
  list("A" = c(120, 140), "B" = c(6, 5), "C" = c(4, 3))
)
ggplot(simdf, aes(time, y, group = interaction(group, id))) +
  geom_line(aes(color = group)) +
  labs(title = "gumbel")

simdf <- growthSim("double logistic",
  n = 20, t = 70,
  params = list(
    "A" = c(200, 160), "B" = c(13, 11), "C" = c(3, 3.5),
    "A2" = c(400, 300), "B2" = c(35, 40), "C2" = c(3.25, 2.75)
  )
)
ggplot(simdf, aes(time, y, group = interaction(group, id))) +
  geom_line(aes(color = group)) +
  labs(title = "Double Logistic")

simdf <- growthSim("double gompertz",
  n = 20, t = 100,
  params = list(
    "A" = c(180, 140), "B" = c(13, 11), "C" = c(0.2, 0.2),
    "A2" = c(400, 300), "B2" = c(50, 50), "C2" = c(0.1, 0.1)
  )
)
ggplot(simdf, aes(time, y, group = interaction(group, id))) +
  geom_line(aes(color = group)) +
  labs(title = "Double Gompertz")

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

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

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

simdf <- growthSim("logarithmic",
  n = 20, t = 25,
  params = list("A" = c(2, 1.7))
)
ggplot(simdf, aes(time, y, group = interaction(group, id))) +
  geom_line(aes(color = group)) +
  labs(title = "Logarithmic")

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

simdf <- growthSim("bragg",
  n = 20, t = 100,
  list("A" = c(10, 15), "B" = c(0.01, 0.02), "C" = c(50, 60))
)
ggplot(simdf, aes(time, y, group = interaction(group, id))) +
  geom_line(aes(color = group)) +
  labs(title = "bragg")

# simulating models from segmented growth models

simdf <- growthSim(
  model = "linear + linear", n = 20, t = 25,
  params = list("linear1A" = c(16, 11), "linear2A" = c(0.75, 0.7), "changePoint1" = c(11, 14))
)
ggplot(simdf, aes(time, y, group = interaction(group, id))) +
  geom_line(aes(color = group)) +
  labs(title = "linear + linear")

simdf <- growthSim(
  model = "linear + linear decay", n = 20, t = 25,
  params = list("linear1A" = c(16, 11), "linear2A" = c(3, 2), "changePoint1" = c(11, 14))
)
ggplot(simdf, aes(time, y, group = interaction(group, id))) +
  geom_line(aes(color = group)) +
  labs(title = "linear + linear decay")

simdf <- growthSim(
  model = "linear + linear + logistic", n = 20, t = 50,
  params = list(
    "linear1A" = c(16, 11), "linear2A" = c(3, 4), # linear slopes, very intuitive
    "changePoint1" = c(11, 14), "changePoint2" = c(10, 12),
    # changepoint1 is standard, changepoint2 happens relative to changepoint 1
    "logistic3A" = c(200, 210), "logistic3B" = c(20, 25), "logistic3C" = c(3, 3)
  )
)
# similar to changepoint2, the asymptote and inflection point are relative to the starting
# point of the logistic growth component. This is different than the model output
# if you were to fit a curve to this model using `growthSS`.
ggplot(simdf, aes(time, y, group = interaction(group, id))) +
  geom_line(aes(color = group)) +
  labs(title = "linear + linear + logistic")

Description

Output from this should be passed to fitGrowth to fit the specified model.

Usage

growthSS(
  model,
  form,
  sigma = NULL,
  df,
  start = NULL,
  pars = NULL,
  type = "brms",
  tau = 0.5,
  hierarchy = NULL
)

Arguments

Details

Default priors are not provided, but these can serve as starting points for each distribution. You are encouraged to use growthSim to consider what kind of trendlines result from changes to your prior and for interpretation of each parameter. The plotPrior function can be used to do prior predictive checks. You should not looking back and forth at your data trying to match your observed growth exactly with a prior distribution, rather this should be informed by an understanding of the plants you are using and expectations based on previous research. For the "double" models the parameter interpretation is the same as for their non-double counterparts except that there are A and A2, etc. It is strongly recommended to familiarize yourself with the double sigmoid distributions using growthSim before attempting to model one. Additionally, those distributions are intended for use with long delays in an experiment, think stress recovery experiments, not for minor hiccups in plant growth.

• Logistic: list('A' = 130, 'B' = 12, 'C' = 3)

• Gompertz: list('A' = 130, 'B' = 12, 'C' = 1.25)

• Weibull: list('A' = 130, 'B' = 2, 'C' = 2)

• Frechet: list('A' = 130, 'B' = 5, 'C' = 6)

• Gumbel: list('A' = 130, 'B' = 6, 'C' = 4)

• Double Logistic: list('A' = 130, 'B' = 12, 'C' = 3, 'A2' = 200, 'B2' = 25, 'C2' = 1)

• Double Gompertz: list('A' = 130, 'B' = 12, 'C' = 0.25, 'A2' = 220, 'B2' = 30, 'C2' = 0.1)

• Monomolecular: list('A' = 130, 'B' = 2)

• Exponential: list('A' = 15, 'B' = 0.1)

• Linear: list('A' = 1)

• Power Law: list('A' = 13, 'B' = 2)

See details below about parameterization for each model option.

• Logistic: 'A / (1 + exp( (B-x)/C) )' Where A is the asymptote, B is the inflection point, C is the growth rate.

• Gompertz: 'A * exp(-B * exp(-C*x))' Where A is the asymptote, B is the inflection point, C is the growth rate.

• Weibull: 'A * (1-exp(-(x/C)^B))' Where A is the asymptote, B is the weibull shape parameter, C is the weibull scale parameter.

• Frechet: 'A * exp(-((x-0)/C)^(-B))' Where A is the asymptote, B is the frechet shape parameter, C is the frechet scale parameter. Note that the location parameter (conventionally m) is 0 in these models for simplicity but is still included in the formula.

• Gumbel: 'A * exp(-exp(-(x-B)/C))' Where A is the asymptote, B is the inflection point (location), C is the growth rate (scale).

• Double Logistic: 'A / (1+exp((B-x)/C)) + ((A2-A) /(1+exp((B2-x)/C2)))' Where A is the asymptote, B is the inflection point, C is the growth rate, A2 is the second asymptote, B2 is the second inflection point, and C2 is the second growth rate.

• Double Gompertz: 'A * exp(-B * exp(-C*x)) + ((A2-A) * exp(-B2 * exp(-C2*(x-B))))' Where A is the asymptote, B is the inflection point, C is the growth rate, A2 is the second asymptote, B2 is the second inflection point, and C2 is the second growth rate.

• Monomolecular: 'A-A * exp(-B * x)' Where A is the asymptote and B is the growth rate.

• Exponential: 'A * exp(B * x)' Where A is the scale parameter and B is the growth rate.

• Linear: 'A * x' Where A is the growth rate.

• Power Law: 'A * x^(B)' Where A is the scale parameter and B is the growth rate.

• Bragg: 'A * exp(-B * (x - C) ^ 2)' This models minima and maxima as a dose-response curve where A is the max response, B is the "precision" or slope at inflection, and C is the x position of the max response.

• Lorentz: 'A / (1 + B * (x - C) ^ 2)' This models minima and maxima as a dose-response curve where A is the max response, B is the "precision" or slope at inflection, and C is the x position of the max response. Generally Bragg is preferred to Lorentz for dose-response curves.

• Beta: 'A * (((x - D) / (C - D)) * ((E - x) / (E - C)) ^ ((E - C) / (C - D))) ^ B' This models minima and maxima as a dose-response curve where A is the Maximum Value, B is a shape/concavity exponent similar to the sum of alpha and beta in a Beta distribution, C is the position of maximum value, D is the minimum position where distribution > 0, E is the maximum position where distribution > 0. This is a difficult model to fit but can model non-symmetric dose-response relationships which may sometimes be worth the extra effort.

Note that for these distributions parameters do not exist in a vacuum. Changing one will make the others look different in the resulting data. The growthSim function can be helpful in familiarizing further with these distributions.

Using the brms backend the sigma argument optionally specifies a sub model to account for heteroskedasticity. Any combination of models (except for decay models) can be specified in the sigma term. If you need variance to raise and lower then a gam/spline is the most appropriate option.

Using the brms backend a model with lots of parameters may be difficult to estimate if there are lots of groups. If you have very many levels of your "group" variable in a complex model then consider fitting models to subsets of the "group" variable and using combineDraws to make a data.frame for hypothesis testing.

Limits on the Y variable can be specified in the brms backend. This should generally be unnecessary and will make the model slower to fit and potentially more difficult to set priors on. If you do have a limited phenotype (besides the normal positive constraint for growth models) then this may be helpful, one situation may be canopy coverage percentage which is naturally bounded at an upper and lower limit. To specify these limits add square brackets to the Y term with upper and lower limits such as "y[0,100] ~ time|id/group". Other "Additional response information" such as resp_weights or standard errors can be specified using the brms backend, with those options documented fully in the brms::brmsformula details.

There are also three supported submodel options for nlme models, but a varFunc object can also be supplied, see ?nlme::varClasses.

• none: varIdent(1|group), which models a constant variance separately for each group.

• power: varPower(x|group), which models variance as a power of x per group.

• exp: varExp(x|group), which models variance as an exponent of x per group.

Survival models can be fit using the "survival" keyword in the model specification. Using the "brms" backend (type argument) you can specify "weibull" (the default) or "binomial" for the distribution to use in that model so that the final model string would be "survival binomial" or "survival weibull" which is equivalent to "survival". Time to event data is very different than standard phenotype data, so the formula argument should include a cutoff for the Y variable to count as an "event". For example, if you were checking germination using area and wanted to use 50 pixels as a germinated plant your formula would be area > 50 ~ time|id/group. Internally the input dataframe will be converted to time-to-event data based on that formula. Alternatively you can make your own time to event data and supply that to growthSS. In that case your data should have columns called "n_events" (number of individuals experiencing the event at this time) and "n_eligible" (number of individuals who had not experienced the event at least up to this time) for the binomial model family OR "event" (binary 1,0 for TRUE, FALSE) for the Weibull model family. Note that since these are linear models using different model families the priors are handled differently. For survival models the default priors are weak regularizing priors (Normal(0,5)) on all parameters. If you wish to specify your own priors you can supply them as brmsprior objects or as a list such as priors = list("group1" = c(0,3), "group2" = c(0,1)) where the order of values is Mu, Sigma. Any non-brms backend will instead use survival::survreg to fit the model unless the "flexsurv" type is specified. Distributions will be passed to survreg where options are "weibull", "exponential", "gaussian", "logistic","lognormal" and "loglogistic" if type = "survreg" or to flexsurv::flexsurvreg if type = "flexsurv" where options are "gengamma", "gengamma.orig", "genf", "genf.orig", "weibull", "gamma", "exp", "llogis", "lnorm", "gompertz", "exponential", and "lognormal". In flexsurvreg distributional modeling is supported and additional formula can be passed as a list to the sigma argument of growthSS in the same way as to the anc argument of flexsurv::flexsurvreg. Further additional arguments should be supplied via fitGrowth if desired.

Value

A named list of elements to make it easier to fit non linear growth models with several R packages.

For brms models the output contains:

formula: A brms::bf formula specifying the growth model, autocorrelation, variance submodel, and models for each variable in the growth model. prior: A brmsprior/data.frame object. initfun: A function to randomly initialize chains using a random draw from a gamma distribution (confines initial values to positive and makes correct number of initial values for chains and groups). df The data input, with dummy variables added if needed and a column to link groups to their factor levels. family The model family, currently this will always be "student". pcvrForm The form argument unchanged. This is returned so that it can be used later on in model visualization. Often it may be a good idea to save the output of this function with the fit model, so having this can be useful later on.

For quantreg::nlrq models the output contains:

formula: An nls style formula specifying the growth model with groups if specified. taus: The quantiles to be fit start: The starting values, typically these will be generated from the growth model and your data in a similar way as shown in stats::selfStart models. df The input data for the model. pcvrForm The form argument unchanged.

For nls models the output is the same as for quantreg::nlrq models but without taus returned.

For nlme::nlme models the output contains:

formula: An list of nlme style formulas specifying the model, fixed and random effects, random effect grouping, and variance model (weights). start: The starting values, typically these will be generated from the growth model and your data in a similar way as shown in stats::selfStart models. df The input data for the model. pcvrForm The form argument unchanged.

For all models the type and model are also returned for simplicity downstream.

See Also

fitGrowth for fitting the model specified by this list and mvSS for the multi-value trait equivalent.

Examples

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", df = simdf,
  start = list("A" = 130, "B" = 12, "C" = 3), type = "brms"
)
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.

fit_test <- fitGrowth(ss,
  iter = 500, cores = 1, chains = 1, backend = "cmdstanr",
  control = list(adapt_delta = 0.999, max_treedepth = 20)
)



# formulas and priors will look different if there is only one group in the data

ex <- growthSim("linear", n = 20, t = 25, params = list("A" = 2))
ex_ss <- growthSS(
  model = "linear", form = y ~ time | id / group, sigma = "spline",
  df = ex, start = list("A" = 1), type = "brms"
)

ex_ss$prior # no coef level grouping for priors
ex_ss$formula # intercept only model for A

ex2 <- growthSim("linear", n = 20, t = 25, params = list("A" = c(2, 2.5)))
ex2_ss <- growthSS(
  model = "linear", form = y ~ time | id / group, sigma = "spline",
  df = ex2, start = list("A" = 1), type = "brms"
)
ex2_ss$prior # has coef level grouping for priors
ex2_ss$formula # specifies an A intercept for each group and splines by group for sigma

Description

EMD can get very heavy with large datasets. For an example lemnatech dataset filtering for images from every 5th day there are 6332^2 = 40,094,224 pairwise EMD values. In long format that's a 40 million row dataframe, which is unwieldy. This function is to help reduce the size of datasets before comparing histograms and moving on with matrix methods or network analysis.

Usage

mv_ag(
  df,
  group,
  mvCols = "frequencies",
  n_per_group = 1,
  outRows = NULL,
  keep = NULL,
  parallel = getOption("mc.cores", 1),
  traitCol = "trait",
  labelCol = "label",
  valueCol = "value",
  id = "image"
)

Arguments

Value

Returns a dataframe summarized by the specified groups over the multi-value traits.

Examples

s1 <- mvSim(
  dists = list(runif = list(min = 15, max = 150)),
  n_samples = 10,
  counts = 1000,
  min_bin = 1,
  max_bin = 180,
  wide = TRUE
)
mv_ag(s1, group = "group", mvCols = "sim_", n_per_group = 2)

Description

mvSim can be used to simulate data for example models/plots.

Usage

mvSim(
  dists = list(rnorm = list(mean = 100, sd = 15)),
  n_samples = 10,
  counts = 1000,
  min_bin = 1,
  max_bin = 180,
  wide = TRUE,
  binwidth = 1,
  t = NULL,
  model = "linear",
  params = list(A = 10)
)

Arguments

Value

Returns a dataframe of example multi-value trait data simulated from specified distributions.

Examples

library(extraDistr) # for rmixnorm
library(ggplot2)
dists <- list(
  rmixnorm = list(mean = c(70, 150), sd = c(15, 5), alpha = c(0.3, 0.7)),
  rnorm = list(mean = 90, sd = 3)
)
x <- mvSim(dists = dists, wide = FALSE)
dim(x)
x2 <- mvSim(dists = dists)
dim(x2)

ggplot(x, aes(
  x = as.numeric(sub("sim_", "", variable)),
  y = value, group = interaction(group, id), fill = group
)) +
  geom_col(position = "identity", alpha = 0.25) +
  pcv_theme() +
  labs(x = "bin")
dists = list(rnorm = list(mean = 30, sd = 15), rnorm = list(mean = 25, sd = 10))
x3 <- mvSim(
  dists = dists, wide = FALSE, # here we make longitudinal data
  t = 10, model = "linear", params = list("A" = c(10, 5))
)
ggplot(x3, aes(
  x = as.numeric(sub("sim_", "", variable)),
  y = value, group = interaction(group, id), fill = group
)) +
  facet_wrap(~times) +
  geom_col(position = "identity", alpha = 0.25) +
  pcv_theme() +
  labs(x = "bin")

Description

This function provides a simplified interface to modeling multi-value traits using growthSS. Output from this should be passed to fitGrowth to fit the specified model.

Usage

mvSS(
  model = "linear",
  form,
  sigma = NULL,
  df,
  start = NULL,
  pars = NULL,
  type = "brms",
  tau = 0.5,
  hierarchy = NULL,
  spectral_index = c("none", "ari", "ci_rededge", "cri550", "cri700", "egi", "evi",
    "gdvi", "mari", "mcari", "mtci", "ndre", "ndvi", "pri", "psnd_chlorophyll_a",
    "psnd_chlorophyll_b", "psnd_caroteniods", "psri", "pssr_chlorophyll_a",
    "pssr_chlorophyll_b", "pssr_caroteniods", "rgri", "rvsi", "savi", "sipi", "sr",
    "vari", "vi_green", "wi", "fvfm", "fqfm")
)

Arguments

Value

A named list of plots showing prior distributions that growthSS would use, optionally with a plot of simulated growth curves using draws from those priors.

See Also

fitGrowth for fitting the model specified by this list.

Examples

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 = "none"
)

mod1 <- fitGrowth(ss1, backend = "cmdstanr", iter = 1000, chains = 1, cores = 1)
growthPlot(mod1, ss1$pcvrForm, df = ss1$df)


# when the model is longitudinal the same model is possible with growthSS

m1 <- mvSim(
  dists = list(
    rnorm = list(mean = 100, sd = 30),
    rnorm = list(mean = 110, sd = 25),
    rnorm = list(mean = 120, sd = 20),
    rnorm = list(mean = 135, sd = 15)
  ),
  wide = FALSE, n = 6
)
m1$time <- rep(1:4, times = 6 * 180)
m2 <- mvSim(
  dists = list(
    rnorm = list(mean = 85, sd = 25),
    rnorm = list(mean = 95, sd = 20),
    rnorm = list(mean = 105, sd = 15),
    rnorm = list(mean = 110, sd = 15)
  ),
  wide = FALSE, n = 6
)
m2$time <- rep(1:4, times = 6 * 180)
mv_df2 <- rbind(m1, m2)
mv_df2$group <- rep(c("a", "b"), each = 4320)
mv_df2 <- mv_df2[mv_df2$value > 0, ]
mv_df2$label <- as.numeric(gsub("sim_", "", mv_df2$variable))
ss_mv0 <- mvSS(
  model = "linear", form = label | value ~ group, df = mv_df2,
  start = list("A" = 50), type = "brms", spectral_index = "ci_rededge"
)
ss_mv0 # non longitudinal model setup

ss_mv1 <- mvSS(
  model = "linear", form = label | value ~ time | group, df = mv_df2,
  start = list("A" = 50), type = "brms", spectral_index = "ci_rededge"
)
ss_mv1
ss_mv2 <- growthSS(
  model = "skew_normal: linear",
  form = label | resp_weights(value) + trunc(lb = -1, ub = Inf) ~ time | group,
  df = mv_df2, start = list("A" = 50)
)
ss_mv2
# ignoring environments and other such details these are identical except for the
# function call.
unlist(lapply(names(ss_mv1), function(nm) {
  if (!identical(ss_mv1[[nm]], ss_mv2[[nm]],
    ignore.environment = TRUE,
    ignore.srcref = TRUE
  )) {
    if (!identical(as.character(ss_mv1[[nm]]), as.character(ss_mv2[[nm]]))) {
      nm
    }
  }
}))


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)
}