Package {biometryassist}

biometryassist: Functions to Assist Design and Analysis of Agronomic Experiments

Description

Provides functions to aid in the design and analysis of agronomic and agricultural experiments through easy access to documentation and helper functions, especially for users who are learning these concepts. While not required for most functionality, this package enhances the ‘asreml' package which provides a computationally efficient algorithm for fitting mixed models using Residual Maximum Likelihood. It is a commercial package that can be purchased as ’ASReml-R' from 'VSNi' https://vsni.co.uk/, who will supply a zip file for local installation/updating (see https://asreml.kb.vsni.co.uk/).

Author(s)

Maintainer: Sam Rogers biometrytraining@adelaide.edu.au

Authors:

• Sam Rogers biometrytraining@adelaide.edu.au

• Sharon Nielsen biometrytraining@adelaide.edu.au

• Annie Conway biometrytraining@adelaide.edu.au

Other contributors:

• Michael Mumford [contributor]

• University of Adelaide (https://adelaide.edu.au/) [copyright holder, funder]

• Grains Research and Development Corporation (https://grdc.com.au/) [copyright holder, funder]

See Also

Useful links:

• https://biometryhub.github.io/biometryassist/

• Report bugs at https://github.com/biometryhub/biometryassist/issues

Add buffers to an existing design

Description

Add buffers to an existing design

Usage

add_buffers(design_obj, type, by = NULL)

Arguments

Value

The modified design object with buffers added

Examples

# Create a simple CRD design
des <- design(type = "crd", treatments = c("A", "B"), reps = 3, nrows = 2, ncols = 3, seed = 42)

# Plot the original design
autoplot(des)

# Add edge buffers (a border around the whole field) and plot
autoplot(add_buffers(des, type = "edge"))

# A buffer row between every pair of rows
autoplot(add_buffers(des, type = "row"))

# Two buffer rows between every pair of rows
autoplot(add_buffers(des, type = "double row"))

# For block designs, buffer the boundaries between blocks
des_rcbd <- design(type = "rcbd", treatments = LETTERS[1:4], reps = 4,
                   nrows = 4, ncols = 4, brows = 2, bcols = 2, seed = 42)

# One buffer at each internal block boundary
autoplot(add_buffers(des_rcbd, type = "block"))

# A full buffer ring surrounding each block
autoplot(add_buffers(des_rcbd, type = "double block"))

# Use `by` to buffer where a grouping column changes. In a split-plot the
# `plots` column identifies each wholeplot, so this rings every wholeplot:
des_sp <- design(type = "split", treatments = c("A", "B"), sub_treatments = 1:4,
                 reps = 4, nrows = 8, ncols = 4, brows = 4, bcols = 2, seed = 42)
autoplot(add_buffers(des_sp, type = "double block", by = "plots"))

# ... or just a buffer row between the wholeplot bands
autoplot(add_buffers(des_sp, type = "row", by = "plots"))

Format the aliased-levels note for the comparison print methods

Description

Shared by print.mct(), print.pairwise_comparisons() and print.reference_comparisons() so all three report aliasing identically. Returns NULL when nothing was aliased. The levels are listed when there are few; once there are more than 6 the list is collapsed to a count (to avoid a large block) and the user is shown how to retrieve them from the aliased attribute.

Usage

aliased_note(aliased)

Arguments

Value

A single-line character string, or NULL.

Generate automatic plots for objects generated in biometryassist

Description

ggplot2::autoplot() methods are provided for the objects created by biometryassist. See the per-class methods for the available options: autoplot.mct() for multiple_comparisons() output, autoplot.design() for design() output, and autoplot.pairwise_comparisons() / autoplot.reference_comparisons(), which are documented alongside their respective functions.

Usage

autoplot(object, ...)

Arguments

Value

A ggplot2 object.

See Also

autoplot.mct(), autoplot.design(), multiple_comparisons() and design()

Plot the layout of an experimental design

Description

Produces a plot of the plot/field layout for a design() result, with plots coloured by treatment and block boundaries drawn for blocked designs.

Usage

## S3 method for class 'design'
autoplot(
  object,
  rotation = 0,
  size = 4,
  margin = FALSE,
  palette = "default",
  row = NULL,
  column = NULL,
  block = NULL,
  treatments = NULL,
  legend = TRUE,
  ...
)

Arguments

Value

A ggplot2 object.

See Also

design()

Examples

des.out <- design(type = "crd", treatments = c(1, 5, 10, 20),
                  reps = 5, nrows = 4, ncols = 5, seed = 42, plot = FALSE)
autoplot(des.out)

# Colour blind friendly colours
autoplot(des.out, palette = "colour-blind")

# Alternative colour scheme
autoplot(des.out, palette = "plasma")

# Custom colour palette
autoplot(des.out, palette = c("#ef746a", "#3fbfc5", "#81ae00", "#c37cff"))

# Visualise different components of a split plot design
des.out <- design(type = "split", treatments = c("A", "B"), sub_treatments = 1:4,
reps = 4, nrows = 8, ncols = 4, brows = 4, bcols = 2, seed = 42)

# Show the wholeplot components
autoplot(des.out, treatments = wholeplots)

# Display block level
autoplot(des.out, treatments = block)

Plot the predicted means from a multiple comparisons test

Description

Produces a plot of the predicted means from a multiple_comparisons() result, with error bars (or a single Tukey's HSD reference bar) and significance-group lettering.

Usage

## S3 method for class 'mct'
autoplot(
  object,
  size = 4,
  label_height = 0.1,
  rotation = 0,
  axis_rotation = rotation,
  label_rotation = rotation,
  type = "point",
  errorbar_type = "ci",
  include_errorbar = TRUE,
  include_lettering = TRUE,
  trans_scale = FALSE,
  ...
)

Arguments

Value

A ggplot2 object.

See Also

multiple_comparisons()

Examples

dat.aov <- aov(Petal.Width ~ Species, data = iris)
output <- multiple_comparisons(dat.aov, classify = "Species")
autoplot(output, label_height = 0.5)

# Join the means with a line
autoplot(output, type = "line", label_height = 0.5)

# Show a single Tukey's HSD reference bar instead of per-mean intervals
autoplot(output, errorbar_type = "hsd")

Pairwise comparisons of predicted means

Description

Test a chosen set of pairwise differences between the predicted means of a fitted model, with multiplicity adjustment over the chosen set and optional splitting into independent subgroups. Unlike multiple_comparisons(), which is means-centric and summarises all pairwise comparisons via a compact letter display, pairwise_comparisons() is difference-centric: it returns a tidy table with one row per requested comparison. This honestly represents selective and irregular comparison sets, for which letter groupings are not valid.

Usage

## S3 method for class 'pairwise_comparisons'
autoplot(object, ..., axis_rotation = 0, label_rotation = 0)

pairwise_comparisons(
  model.obj,
  classify,
  pairs = NULL,
  contrasts = NULL,
  adjust = "holm",
  by = NULL,
  sig = 0.05,
  include_means = TRUE,
  descending = NULL,
  ...
)

## S3 method for class 'pairwise_comparisons'
print(x, decimals = 2, ...)

Arguments

Details

Relationship to the other comparison functions

pairwise_comparisons(), multiple_comparisons() and reference_comparisons() share the same predicted means and standard errors of differences; they differ in which comparisons they report. Testing all pairs here (pairs = NULL) with a given adjust yields the same adjusted p-values as multiple_comparisons() with that adjust, while reference_comparisons() is the special case of comparing every level against a single control. Tukey's HSD is exact only for the complete set of all pairwise comparisons, so adjust = "tukey" is not valid for pairwise_comparisons() (use multiple_comparisons()).

pairs syntax and sign convention

The estimate for a pair is mean(level1) - mean(level2), in the order written ("A-B" gives A - B). The level1 and level2 columns make the sign unambiguous. With : joining interaction cells and - separating the two sides of a pair, a level name may itself contain -: each - is tried as the separator and the split that yields two valid levels is used (so "A-D-xyz" resolves to A versus D-xyz when those are the real levels). Only when more than one split is valid is the pair genuinely ambiguous, and the list form is then required (the error shows the candidate splits). Reversed or duplicated pairs ("A-B" and "B-A") are de-duplicated with a warning, since duplicates would inflate the adjustment family.

For an interaction classify, the :-joined components of each level label must be in the same order as classify (e.g. "A:X" for classify = "Trt:Site", not "X:A"); the same applies to the level names used in contrasts. A label that names no existing cell is rejected with the list of available levels, but a mis-ordered label that happens to name another valid cell is used silently — so when the factors share level names, check the component order matches classify.

by semantics

by must be a subset of the classify factors. Within each group, pair labels reference the remaining (non-by) factor levels. For example, classify = "Trt:Site", by = "Site", pairs = "A-B" compares Trt levels A and B within each Site. A group with fewer than two levels is skipped with a warning. In an unbalanced design a requested comparison may reference a level that is absent from some groups: such a comparison is skipped (with a warning) in the groups where it cannot be computed and reported in those where it can. A level that is absent from every group — or that was aliased — is instead an error, since that indicates a mistake rather than an incomplete design.

Standard error and degrees of freedom for contrasts

The contrast variance is ⁠c' V c⁠, where V is the variance-covariance matrix of the predicted means taken directly from the fitting engine (no reconstruction). The degrees of freedom depend on the engine:

• For models predicted via emmeans (aov, lm, nlme::lme(), lme4::lmer(), and aov with Error() strata) the estimate, standard error and degrees of freedom are obtained directly from emmeans::contrast() on the model's reference grid. The degrees of freedom are therefore the exact contrast df for that engine (Satterthwaite or Kenward-Roger for mixed models, containment for aov/Error()), including for contrasts spanning more than two levels.

• For asreml models V is the prediction error covariance from predict(..., vcov = TRUE), and the degrees of freedom are the term's denominator df from asreml::wald(denDF = "default") (a single Kenward-Roger-style value shared by all contrasts within the term, as ASReml-R does not provide a per-contrast approximate df).

Transformations

Comparisons are reported on the model scale. A difference of transformed means does not back-transform to a difference on the original scale (for log/logit it is a ratio), so back-transformation is not attempted; a warning is issued if the response appears to be transformed in the model formula.

Confidence intervals

The per-comparison confidence intervals are not simultaneity-adjusted: as with the confidence-interval/letter note in multiple_comparisons(), an interval may exclude zero while the adjusted p-value is ⁠>= sig⁠. When this happens a one-time message() notes it (suppressible with suppressMessages()).

Value

autoplot.pairwise_comparisons() returns a ggplot2 object: a forest plot of the estimated differences with their confidence intervals and a dashed reference line at zero, faceted by the by variable(s) when present. Comparisons that are significant at the adjusted sig level are flagged with an asterisk (*) — prefixed to the y-axis label when unfaceted (keeping the labels right-justified against the axis), or beside the interval when faceted by by.

A data.frame of class pairwise_comparisons with one row per (group × comparison) and columns: any by column(s), level1, level2, comparison, estimate, (optionally level1.mean and level2.mean), std.error, statistic, df, p.value (adjusted), conf.low and conf.high. Stored at full precision; rounding for display is controlled by print.pairwise_comparisons(). For the contrasts form, the level1/level2 and mean columns are omitted and comparison holds the contrast name. If any levels were aliased (not estimable) in the model they are dropped from the comparisons (with a warning) and recorded in an aliased attribute.

print.pairwise_comparisons() invisibly returns x.

Supported model types

The comparison functions (multiple_comparisons(), pairwise_comparisons() and reference_comparisons()) work with any model for which a get_predictions() method is defined. These are currently:

ARTool (art) models are supported by resplot() but not by the comparison functions: the aligned rank transform makes mean-based comparisons inappropriate. Use ARTool::art.con() for contrasts on ART models instead.

sommer mmer models (the legacy interface) are supported by resplot() but not by the comparison functions: current sommer provides no predict() method for mmer. Refit with sommer::mmes() to use the comparison functions.

To add a new engine, write a ⁠get_predictions.<class>()⁠ method returning a list with elements predictions, sed, df, ylab and aliased_names (plus emmeans_grid for engines backed by emmeans::emmeans()), and add a row to the table above.

See Also

multiple_comparisons() for means-and-letters output, and reference_comparisons() for comparing every level against a single control. For guidance on choosing between them and on multiplicity adjustments, see vignette("choosing-multiple-comparisons", "biometryassist").

Examples

# Forest plot of pairwise differences (significant comparisons marked with *)
dat.aov <- aov(Petal.Width ~ Species, data = iris)
pc <- pairwise_comparisons(dat.aov, classify = "Species")
autoplot(pc)
dat.aov <- aov(Petal.Width ~ Species, data = iris)

# All pairwise comparisons (Holm-adjusted by default)
pairwise_comparisons(dat.aov, classify = "Species")

# A selected subset
pairwise_comparisons(
    dat.aov,
    classify = "Species",
    pairs = c("setosa-versicolor", "setosa-virginica")
)

# A general (non-pairwise) contrast: setosa vs the average of the others
pairwise_comparisons(
    dat.aov,
    classify = "Species",
    contrasts = list(
        "setosa vs rest" = c(setosa = 1, versicolor = -0.5, virginica = -0.5)
    )
)

# `by`: the same comparisons, adjusted independently within each group. Here a
# 2 x 3 factorial - compare tension levels within each wool type. The `pairs`
# labels reference the remaining (tension) factor.
m_wb <- aov(breaks ~ wool * tension, data = warpbreaks)
pairwise_comparisons(
    m_wb,
    classify = "wool:tension",
    by = "wool",
    pairs = c("L-M", "L-H")
)

Compare predicted means against a single reference (control) level

Description

Compare every level of a treatment factor against one chosen reference (control) level, returning a tidy, means-centric table: the predicted mean of each level, the reference mean, their difference, and a multiplicity-adjusted p-value for the difference. This is the honest representation of the common "how does each treatment compare to the control?" question (a Dunnett-style analysis), where significance attaches cleanly to each treatment because every comparison shares the same reference.

Usage

## S3 method for class 'reference_comparisons'
autoplot(object, ..., axis_rotation = 0, label_rotation = 0)

reference_comparisons(
  model.obj,
  classify,
  reference,
  adjust = "dunnett",
  by = NULL,
  sig = 0.05,
  include_means = TRUE,
  descending = NULL,
  ...
)

## S3 method for class 'reference_comparisons'
print(x, decimals = 2, ...)

Arguments

Details

Unlike multiple_comparisons() (all pairs, means + letters) and pairwise_comparisons() (selected differences), reference_comparisons() compares each level to a single control and presents the result around the means. It works for every model supported by multiple_comparisons().

Why Dunnett is the default

Comparing several treatments to one control has an exact simultaneous procedure — Dunnett's test — which uses the joint multivariate-t distribution of the comparisons (accounting for the correlation induced by the shared control). It is more powerful than applying a generic adjustment, so it is the default. The exact test requires a single (common) degrees-of-freedom; for models that report comparison-specific df, reference_comparisons() falls back to "holm" with a warning. Any stats::p.adjust.methods method may be requested explicitly (a message notes that Dunnett is the exact option).

The multivariate-t routine (mvtnorm::pmvt()) requires an integer degrees-of-freedom, so a fractional denominator df (such as an ASReml-R Kenward-Roger denDF) is rounded to the nearest integer for the Dunnett calculation when two or more comparisons are made. The reported df column is the exact (unrounded) value, and a single comparison uses the exact df. In practice this rounding only affects asreml models with a fractional denDF (aov/lm/lme have integer df, and mixed models with comparison-specific df fall back to Holm).

The Dunnett correlation structure is taken from the variance-covariance matrix of the predicted means that the prediction machinery returns, so the exact test is available for every supported model engine. With adjust = "dunnett" the confidence intervals are the simultaneous Dunnett intervals and therefore agree with the adjusted test (an interval excludes zero exactly when the comparison is significant). With a stats::p.adjust() method the intervals are per-comparison and may disagree with the adjusted p-value.

by semantics

by must be a subset of the classify factors. Within each group, the reference and the compared levels reference the remaining (non-by) factor. For example classify = "Trt:Site", by = "Site", reference = "Control" compares every Trt level against Control within each Site, adjusted within each Site. A group missing the reference, or with fewer than two levels, is skipped with a warning.

To compare a control level against the others within every combination of the remaining factors, put all the other factors in by. The within-group factor is then the single control factor, so reference is just that level (e.g. "0"). For example, with classify = "time:variety:dose" and a zero dose labelled "0", by = c("time", "variety"), reference = "0" compares each dose against the zero dose within every time-by-variety cell. Specified this way reference is a single level, so it is unaffected by the order of the factors in classify or by (those orderings only change the column order and group labels of the output, not the comparisons).

Without by, by contrast, reference is the full :-joined cell label and its components must follow the classify order. A mis-ordered label that happens to name another valid cell is used silently (no error), so it is worth double-checking the order matches; the order-proof by form above avoids this entirely.

Transformations

Comparisons are reported on the model scale (a difference of transformed means does not back-transform to a difference on the original scale); a warning is issued if the response appears to be transformed in the model formula.

Value

autoplot.reference_comparisons() returns a ggplot2 object: a means plot with one point per level at its predicted mean, a dashed reference line at the reference mean (marked with a diamond), and an interval around each mean showing the (adjusted) confidence interval for the difference from the reference — so the interval clears the reference line exactly when the comparison is significant (with adjust = "dunnett"). Faceted by the by variable(s) when present. Significant comparisons are flagged with an asterisk (*), prefixed to the y-axis label when unfaceted or beside the interval when faceted.

A data.frame of class reference_comparisons with one row per (group × non-reference level) and columns: any by column(s), level1, level2 (the reference), comparison, estimate, level1.mean, level2.mean (the reference mean), std.error, statistic, df, p.value (adjusted), conf.low and conf.high. The reference level is not given its own row (its mean appears as level2.mean on every row, and in the reference attribute). Stored at full precision; rounding for display is controlled by print.reference_comparisons(). If any levels were aliased (not estimable) in the model they are dropped (with a warning) and recorded in an aliased attribute; an aliased reference is an error.

print.reference_comparisons() invisibly returns x.

Supported model types

The comparison functions (multiple_comparisons(), pairwise_comparisons() and reference_comparisons()) work with any model for which a get_predictions() method is defined. These are currently:

ARTool (art) models are supported by resplot() but not by the comparison functions: the aligned rank transform makes mean-based comparisons inappropriate. Use ARTool::art.con() for contrasts on ART models instead.

sommer mmer models (the legacy interface) are supported by resplot() but not by the comparison functions: current sommer provides no predict() method for mmer. Refit with sommer::mmes() to use the comparison functions.

To add a new engine, write a ⁠get_predictions.<class>()⁠ method returning a list with elements predictions, sed, df, ylab and aliased_names (plus emmeans_grid for engines backed by emmeans::emmeans()), and add a row to the table above.

See Also

multiple_comparisons() for all-pairs means and letters, pairwise_comparisons() for selected differences. For guidance on choosing between them and on multiplicity adjustments, see vignette("choosing-multiple-comparisons", "biometryassist").

Examples

# Means plot of each level vs the reference (significant ones marked with *)
dat.aov <- aov(Petal.Width ~ Species, data = iris)
rc <- reference_comparisons(dat.aov, classify = "Species", reference = "setosa")
autoplot(rc)
dat.aov <- aov(weight ~ feed, data = chickwts)

# Compare every feed against the "casein" control (exact Dunnett)
reference_comparisons(dat.aov, classify = "feed", reference = "casein")

# A different adjustment can be requested (a message notes that Dunnett is the
# exact option for this family):
reference_comparisons(
    dat.aov,
    classify = "feed",
    reference = "casein",
    adjust = "holm"
)

# `by`: compare each level against the reference *within* each group, adjusted
# independently. Here a 2 x 3 factorial - each tension level vs the "L"
# reference within each wool type; autoplot() facets by wool.
m_wb <- aov(breaks ~ wool * tension, data = warpbreaks)
rc_by <- reference_comparisons(
    m_wb,
    classify = "wool:tension",
    reference = "L",
    by = "wool"
)
rc_by
autoplot(rc_by)

Deprecated functions in package biometryassist.

Description

The functions listed below are deprecated and will be removed in the future. When possible, alternative functions with similar functionality are also mentioned. Help pages for deprecated functions are available at help("<function>-deprecated").

Usage

resplt(
  model.obj,
  shapiro = TRUE,
  call = FALSE,
  label.size = 10,
  axes.size = 10,
  call.size = 9,
  mod.obj
)

Value

No return value.

A list containing ggplot2 objects which are diagnostic plots. For resplt, use resplot().

resplt

Residual plots of linear models.

Check if classify term exists in model terms

Description

Checks if the classify term exists in the model terms, handling interaction terms specified in any order (e.g., B:A when model has A:B, or A:C:B when model has A:B:C).

Usage

check_classify_in_terms(classify, model_terms)

Arguments

Value

The classify term as it appears in the model (potentially reordered), or throws an error if not found

Combine plots into final output

Description

Combine plots into final output

Usage

combine_plots(plots, shapiro_result, model_call, call, call.size, label.size)

Arguments

Function to compare package version for mocking

Description

Function to compare package version for mocking

Usage

compare_version(a, b)

Arguments

Value

Numeric. 0 if the numbers are equal, -1 if b is later and 1 if a is later

Create buffers for design plots

Description

Create buffers for design plots

Usage

create_buffers(design, type, by = NULL)

Arguments

Value

The original data frame, updated to include buffers

Create diagnostic plots

Description

Create diagnostic plots

Usage

create_diagnostic_plots(group_residuals, axes.size, label.size)

Arguments

Create the folder macOS needs for licensing

Description

ASReml-R uses Reprise licence management which requires the folder /Library/Application Support/Reprise/ to exist on macOS Big Sur (11) and later. This function checks for the folder and attempts to create it if missing, first without elevated privileges, then via a native macOS authentication dialog using AppleScript.

Usage

create_mac_folder(reprise_path = "/Library/Application Support/Reprise/")

Arguments

Value

Logical TRUE if the folder exists or was successfully created. Stops with an informative error if the folder could not be created and the user cancelled the authentication dialog.

Produce a graph of design layout, skeletal ANOVA table and data frame with complete design

Description

Deprecated: des_info() has been superseded by design(). Please use design() instead. This function will be removed in a future version.

Usage

des_info(
  design.obj,
  nrows,
  ncols,
  brows = NA,
  bcols = NA,
  byrow = TRUE,
  fac.names = NULL,
  fac.sep = c("", " "),
  buffer = NULL,
  plot = TRUE,
  rotation = 0,
  size = 4,
  margin = FALSE,
  save = FALSE,
  savename = paste0(design.obj$parameters$design, "_design"),
  plottype = "pdf",
  return.seed = TRUE,
  quiet = FALSE,
  ...
)

Arguments

Value

A list containing a data frame with the complete design, a ggplot object with plot layout, the seed (if return.seed = TRUE), and the satab object, allowing repeat output of the satab table via cat(output$satab).

Create a complete experimental design with graph of design layout and skeletal ANOVA table

Description

Create a complete experimental design with graph of design layout and skeletal ANOVA table

Usage

design(
  type,
  treatments,
  reps,
  nrows,
  ncols,
  brows = NA,
  bcols = NA,
  byrow = TRUE,
  sub_treatments = NULL,
  fac.names = NULL,
  fac.sep = c("", " "),
  buffer = NULL,
  plot_numbers = FALSE,
  plot = TRUE,
  rotation = 0,
  size = 4,
  margin = FALSE,
  save = FALSE,
  savename = paste0(type, "_design"),
  plottype = "pdf",
  seed = TRUE,
  quiet = FALSE,
  ...
)

Arguments

Details

Supported designs are Completely Randomised (crd), Randomised Complete Block (rcbd), Latin Square (lsd), split-plot (split), strip-plot (strip), and crossed factorial designs via ⁠crossed:<base>⁠ where ⁠<base>⁠ is crd, rcbd, or lsd (e.g. crossed:crd).

If save = TRUE (or "both"), both the plot and the workbook will be saved to the current working directory, with filename given by savename. If one of either "plot" or "workbook" is specified, only that output is saved. If save = FALSE (the default, or equivalently "none"), nothing will be output.

fac.names can be supplied to provide more intuitive names for factors and their levels in factorial and split plot designs. They can be specified in a list format, for example fac.names = list(A_names = c("a", "b", "c"), B_names = c("x", "y", "z")). This will result a design output with a column named A_names with levels ⁠a, b, c⁠ and another named B_names with levels ⁠x, y, z⁠. Labels can also be supplied as a character vector (e.g. c("A", "B")) which will result in only the treatment column names being renamed. Only the first two elements of the list will be used, except in the case of a 3-way factorial design.

... allows extra arguments to be passed to ggsave() for output of the plot. The details of possible arguments can be found in ggplot2::ggsave().

Value

A list containing a data frame with the complete design (⁠$design⁠), a ggplot object with plot layout (⁠$plot.des⁠), the seed (⁠$seed⁠, if return.seed = TRUE), and the satab object (⁠$satab⁠), allowing repeat output of the satab table via cat(output$satab).

Examples

# Completely Randomised Design
des.out <- design(type = "crd", treatments = c(1, 5, 10, 20),
                  reps = 5, nrows = 4, ncols = 5, seed = 42)

# Randomised Complete Block Design
des.out <- design("rcbd", treatments = LETTERS[1:11], reps = 4,
                  nrows = 11, ncols = 4, brows = 11, bcols = 1, seed = 42)

# Latin Square Design
# Doesn't require reps argument
des.out <- design(type = "lsd", c("S1", "S2", "S3", "S4"),
                  nrows = 4, ncols = 4, seed = 42)

# Factorial Design (Crossed, Completely Randomised)
des.out <- design(type = "crossed:crd", treatments = c(3, 2),
                  reps = 3, nrows = 6, ncols = 3, seed = 42)

# Factorial Design (Crossed, Completely Randomised), renaming factors
des.out <- design(type = "crossed:crd", treatments = c(3, 2),
                  reps = 3, nrows = 6, ncols = 3, seed = 42,
                  fac.names = list(N = c(50, 100, 150),
                                   Water = c("Irrigated", "Rain-fed")))

# Factorial Design (Crossed, Randomised Complete Block Design),
# changing separation between factors
des.out <- design(type = "crossed:rcbd", treatments = c(3, 2),
                  reps = 3, nrows = 6, ncols = 3,
                  brows = 6, bcols = 1,
                  seed = 42, fac.sep = c(":", "_"))

# Factorial Design (Nested, Latin Square)
trt <- c("A1", "A2", "A3", "A4", "B1", "B2", "B3")
des.out <- design(type = "lsd", treatments = trt,
                  nrows = 7, ncols = 7, seed = 42)

# Split plot design
des.out <- design(type = "split", treatments = c("A", "B"), sub_treatments = 1:4,
                  reps = 4, nrows = 8, ncols = 4, brows = 4, bcols = 2, seed = 42)

# Alternative arrangement of the same design as above
des.out <- design(type = "split", treatments = c("A", "B"), sub_treatments = 1:4,
                  reps = 4, nrows = 8, ncols = 4, brows = 4, bcols = 2,
                  byrow = FALSE, seed = 42)

# Strip plot design
des.out <- design(type = "strip", treatments = c("A", "B", "C"), sub_treatments = 1:4,
                  reps = 4, nrows = 12, ncols = 4, brows = 3, bcols = 4, seed = 42)

Detect Linux distribution and version

Description

Reads /etc/os-release and extracts distro ID and major version.

Usage

detect_linux()

Value

list(id, like, major)

Detect macOS version

Description

Detect macOS version

Usage

detect_macos()

Download asreml package file

Description

Download asreml package file

Usage

download_asreml_package(url, verbose = FALSE)

Arguments

Export Experimental Design Layout to Excel

Description

Converts an experimental design dataframe into a spatial layout matrix and exports to Excel with optional colour coding by treatment.

Usage

export_design_to_excel(
  design,
  filename = "experimental_design.xlsx",
  value_column = "treatments",
  row = NULL,
  column = NULL,
  palette = "default"
)

Arguments

Details

This function takes an experimental design in long format (with row/col coordinates) and converts it to a matrix layout that matches the spatial arrangement of the experiment, then exports to Excel with formatting and optional colour coding.

Valid palette options include:

• "default" - Spectral palette

• ColorBrewer palettes: "brbg", "piyg", "prgn", "puor", "rdbu", "rdgy", "rdylbu", "rdylgn", "spectral", "set3", "paired"

• Viridis palettes: "viridis", "magma", "inferno", "cividis", "plasma", "rocket", "mako", "turbo"

• colour blind friendly: "colour blind", "color blind", "cb" (uses viridis)

• Custom vector of colours (must match number of unique treatments)

Requires the 'openxlsx2' package to be installed.

Value

Invisibly returns the layout dataframe

Examples

## Not run: 
# Export with default colours
export_design_to_excel(my_design, "treatments", "my_design.xlsx")

# Export without colours
export_design_to_excel(my_design, "treatments", "my_design.xlsx", palette = NULL)

# Export with custom palette
export_design_to_excel(my_design, "treatments", "my_design.xlsx", palette = "viridis")

## End(Not run)

Internal model-information extraction for resplot()

Description

extract_model_info() is the internal generic that resplot() uses to pull the residuals, fitted values and (optionally) the model call from a fitted model. It dispatches on the class of model.obj. It is not exported and is not called directly by users; support for a new model engine is added by writing a new extract_model_info() method.

Usage

extract_model_info(model.obj, call = FALSE)

Arguments

Value

A list with elements facet, facet_name, resids, fits, k and model_call.

Supported model types

resplot() produces residual diagnostics for any model with an extract_model_info() method. These are currently:

This set differs slightly from the comparison functions (see get_predictions()): resplot() additionally supports ARTool (art) models and sommer's legacy mmer interface. Neither is available for the comparison functions — ART uses aligned ranks (use ARTool::art.con()), and current sommer provides no predict() for mmer (refit with sommer::mmes()).

To add a new engine, write an ⁠extract_model_info.<class>()⁠ method returning a list with elements facet, facet_name, resids, fits, k and model_call, and add a row to the table above.

See Also

resplot()

Fetch the ASReml build manifest

Description

Downloads and parses the JSON manifest used to decide which ASReml build (URL + version) matches the current OS and R version.

Usage

fetch_manifest(
  manifest_url =
    "https://raw.githubusercontent.com/biometryhub/biometryassist/main/inst/manifest.json"
)

Arguments

Value

A list as returned by jsonlite::fromJSON() with at least a packages element (a list of package entries).

Find existing asreml file

Description

Find existing asreml file

Usage

find_existing_package()

Find the manifest entry matching the current system

Description

Looks up the manifest for a build exactly matching the slug built by get_r_os(). If no exact match is found, falls back to the highest available OS version that does not exceed the current OS version, for the same OS, R version, and architecture. This handles cases such as running on a newer OS than VSNi currently builds for.

Usage

find_package(manifest, os_ver, warn = TRUE)

Arguments

Details

Windows has no OS version component and uses exact matching only.

Value

A single package entry list from manifest$packages, or NULL if no suitable match was found.

Format final output based on facet structure

Description

Format final output based on facet structure

Usage

format_output_resplot(
  output,
  facet,
  facet_name,
  onepage,
  onepage_cols,
  label.size
)

Arguments

Internal prediction extraction for the comparison functions

Description

get_predictions() is the internal generic that multiple_comparisons(), pairwise_comparisons() and reference_comparisons() use to obtain the predicted means, the standard-error-of-differences (SED) matrix and the degrees of freedom from a fitted model. It dispatches on the class of model.obj. It is not exported and is not called directly by users; support for a new model engine is added by writing a new get_predictions() method.

Usage

get_predictions(model.obj, classify, pred.obj = NULL, ...)

Arguments

Value

A list with elements predictions, sed, df, ylab and aliased_names (and emmeans_grid for emmeans-backed engines).

Supported model types

The comparison functions (multiple_comparisons(), pairwise_comparisons() and reference_comparisons()) work with any model for which a get_predictions() method is defined. These are currently:

ARTool (art) models are supported by resplot() but not by the comparison functions: the aligned rank transform makes mean-based comparisons inappropriate. Use ARTool::art.con() for contrasts on ART models instead.

sommer mmer models (the legacy interface) are supported by resplot() but not by the comparison functions: current sommer provides no predict() method for mmer. Refit with sommer::mmes() to use the comparison functions.

To add a new engine, write a ⁠get_predictions.<class>()⁠ method returning a list with elements predictions, sed, df, ylab and aliased_names (plus emmeans_grid for engines backed by emmeans::emmeans()), and add a row to the table above.

See Also

multiple_comparisons(), pairwise_comparisons(), reference_comparisons()

Detect operating system and construct OS/version key

Description

Determines OS, OS version (if applicable), architecture, and combines this with the R version for downstream download logic.

Usage

get_r_os()

Value

A list with os_ver, os, ver, and arm

Get R major-minor version without dot

Description

Converts R version to a compact form (e.g. 4.4.x -> "44").

Usage

get_r_version_compact()

Value

Character scalar R version

Handle deprecated parameters

Description

Simple internal function to warn about deprecated parameters

Usage

handle_deprecated_param(
  old_param,
  new_param = NULL,
  custom_msg = NULL,
  call_env = parent.frame()
)

Arguments

Value

Nothing, called for side effects (warnings)

Produce a heatmap of variables in a grid layout.

Description

This function plots a heatmap of variables in a grid layout, optionally grouping them.

Usage

heat_map(
  data,
  value,
  x_axis,
  y_axis,
  grouping = NULL,
  raster = TRUE,
  smooth = FALSE,
  palette = "default",
  ...
)

Arguments

Value

A ggplot2 object.

Examples

set.seed(42)
dat <- expand.grid(x = 1:5, y = 1:6)
dat$value <- rnorm(30)
dat$groups <- sample(rep(LETTERS[1:6], times = 5))

heat_map(dat, value, x, y)

# Column names can be quoted, but don't need to be.
heat_map(dat, "value", "x", "y", "groups")

# Different palettes are available
heat_map(dat, value, x, y, palette = "Spectral")


# Arguments in ... are passed through to facet_wrap
heat_map(dat, value, x, y, groups, labeller = ggplot2:::label_both)
heat_map(dat, value, x, y, groups, scales = "free_y")
heat_map(dat, value, x, y, groups, nrow = 1)

Install or update the ASReml-R package

Description

Helper functions for installing or updating the ASReml-R package, intended to reduce the difficulty of finding the correct version for your operating system and R version.

Usage

install_asreml(
  library = .libPaths()[1],
  quiet = FALSE,
  force = FALSE,
  keep_file = FALSE,
  check_version = TRUE
)

update_asreml(...)

Arguments

Details

The ASReml-R package file is downloaded from a shortlink, and if keep_file is TRUE, the package archive file will be saved in the current directory. If a valid path is provided in keep_file, the file will be saved to that path, but all directories are assumed to exist and will not be created. If keep_file does not specify an existing, valid path, an error will be shown after package installation.

Value

Silently returns TRUE if asreml installed successfully or already present, FALSE otherwise. Optionally prints a confirmation message on success.

Examples

## Not run: 
# Example 1: download and install asreml
install_asreml()

# Example 2: install asreml and save file for later
install_asreml(keep_file = TRUE)

# Example 3: install with verbose debugging
install_asreml(quiet = "verbose")

## End(Not run)

Install the ASReml package

Description

Install the ASReml package

Usage

install_asreml_package(save_file, library, quiet, os, verbose = FALSE)

Arguments

Value

TRUE if successful, FALSE otherwise

Install required dependencies

Description

Install required dependencies

Usage

install_dependencies(quiet, library, verbose = FALSE)

Arguments

Convert column number to Excel column letter

Description

Convert column number to Excel column letter

Usage

int2col(num)

Arguments

Value

Character Excel column letter(s)

Detect whether the system is ARM-based

Description

x86_64 is treated as the implicit default.

Usage

is_arm()

Value

Logical indicating ARM architecture

Determine if a Colour is Light

Description

Internal helper function to determine whether a colour is light or dark for appropriate font colour selection (black text on light backgrounds, white text on dark backgrounds).

Usage

is_light_colour(colour)

Arguments

Details

Uses standard luminance calculation: 0.299R + 0.587G + 0.114*B, normalized to 0-1 scale. Coefficients reflect human eye sensitivity to different colours (green > red > blue).

Value

Logical. TRUE if the colour is light (luminance > 0.5), FALSE if dark.

List available biometryassist templates

Description

list_templates() returns a character vector of available analysis templates included with the biometryassist package.

Usage

list_templates()

Value

character() Vector of available template file names.

See Also

use_template() to copy and use a template

Examples

# See what templates are available
list_templates()

Conduct a log-likelihood test for comparing terms in ASReml-R models

Description

Conduct a log-likelihood test for comparing terms in ASReml-R models

Usage

logl_test(
  model.obj,
  rand.terms = NULL,
  resid.terms = NULL,
  decimals = 3,
  numeric = FALSE,
  quiet = FALSE
)

Arguments

Value

A data frame of terms and corresponding log-likelihood ratio test p-values.

Build per-row treatment labels from one or more factor columns

Description

Shared by multiple_comparisons() and pairwise_comparisons() to turn the classify factor column(s) of a predictions data frame into a single label per row. A single factor is returned as-is (preserving its type); multiple factors (an interaction) are joined with sep. The caller chooses the separator (multiple_comparisons() uses "_", pairwise_comparisons() uses ":") and any further processing (e.g. coercion to character).

Usage

make_treatment_labels(pp, vars, sep)

Arguments

Value

A vector of labels, one per row of pp.

Manage the downloaded file

Description

Manage the downloaded file

Usage

manage_file(save_file, keep_file, filename, verbose = FALSE)

Arguments

Value

logical; TRUE if file successfully handled, FALSE otherwise

Map Linux distro to supported download target

Description

Map Linux distro to supported download target

Usage

map_linux_target(info)

Perform Multiple Comparison Tests on a statistical model

Description

A function for comparing and ranking predicted means with Tukey's Honest Significant Difference (HSD) Test.

Usage

multiple_comparisons(
  model.obj,
  classify,
  sig = 0.05,
  int.type = "ci",
  trans = NULL,
  offset = NULL,
  power = NULL,
  decimals = 2,
  descending = FALSE,
  groups = TRUE,
  adjust = "tukey",
  by = NULL,
  plot = FALSE,
  label_height = 0.1,
  rotation = 0,
  save = FALSE,
  savename = "predicted_values",
  ...
)

Arguments

Details

Offset

Some transformations require that data has a small offset to be applied, otherwise it will cause errors (for example taking a log of 0, or the square root of negative values). In order to correctly reverse this offset, if the trans argument is supplied, a value should also be supplied in the offset argument. By default the function assumes no offset was required for a transformation, implying a value of 0 for the offset argument. If an offset value is provided, use the same value as provided in the model, not the inverse. For example, if adding 0.1 to values for a log transformation, add 0.1 in the offset argument.

Power

The power argument allows the specification of arbitrary powers to be back transformed, if they have been used to attempt to improve normality of residuals.

P-value adjustment (adjust)

By default (adjust = "tukey") the function uses Tukey's HSD, which is exact for the complete set of all pairwise comparisons. Alternatively, adjust may be any method accepted by stats::p.adjust(). In that case a matrix of raw two-sided t-test p-values is computed first and the chosen adjustment is applied to the lower-triangle vector of those raw p-values, avoiding the double-adjustment that would occur if Tukey-adjusted p-values were passed to p.adjust(). Bonferroni and Holm control the family-wise error rate and are valid under arbitrary dependence; BH/BY (FDR) are valid under positive dependence, which pairwise comparisons satisfy. The returned ⁠$pairwise_pvalues⁠ always holds the adjusted p-values for the chosen method (the Tukey p-values when adjust = "tukey"). When adjust is not "tukey", ⁠$hsd⁠ is NULL, as an HSD value is only meaningful for Tukey's test.

Grouped comparisons (by)

When by is supplied, the predictions are split into subgroups defined by the level(s) of the by variable(s) and the comparison procedure is run independently within each subgroup. Each subgroup is treated as a separate family of comparisons: there is no pooling and no cross-group p-value adjustment, and letter groupings restart within each subgroup. When by is used, ⁠$pairwise_pvalues⁠ (and ⁠$hsd⁠ where applicable) are returned as named lists with one element per subgroup, and autoplot() facets the plot by the by variable(s) by default. by must leave at least one classify factor to compare within each subgroup, so a single-factor classify cannot be split with by.

Confidence Intervals & Comparison Intervals

The function provides several options for confidence intervals via the int.type argument:

• ci (default): Traditional confidence intervals for individual means. These estimate the precision of each individual mean but may not align with the multiple comparison results. Non-overlapping traditional confidence intervals do not necessarily indicate significant differences in multiple comparison tests.

• tukey: Tukey comparison intervals that are consistent with the multiple comparison test. These intervals are wider than regular confidence intervals and are designed so that non-overlapping intervals correspond to statistically significant differences in the Tukey HSD test. This ensures visual consistency between the intervals and letter groupings.

• ⁠1se⁠ and ⁠2se⁠: Intervals of ±1 or ±2 standard errors around each mean.

• none: No confidence intervals will be calculated or displayed in plots.

By default, the function displays regular confidence intervals (int.type = "ci"), which estimate the precision of individual treatment means. However, when performing multiple comparisons, these regular confidence intervals may not align with the letter groupings from Tukey's HSD test. Specifically, you may observe non-overlapping confidence intervals for treatments that share the same letter group (indicating no significant difference).

This occurs because regular confidence intervals and Tukey's HSD test serve different purposes:

• Regular confidence intervals estimate individual mean precision

• Tukey's HSD controls the family-wise error rate across all pairwise comparisons

To resolve this visual inconsistency, you can use Tukey comparison intervals (int.type = "tukey"). These intervals are specifically designed for multiple comparisons and will be consistent with the letter groupings: non-overlapping Tukey intervals indicate significant differences, while overlapping intervals suggest no significant difference.

The function will issue a message if it detects potential inconsistencies between non-overlapping confidence intervals and letter groupings, suggesting the use of Tukey intervals for clearer interpretation. For multiple comparison contexts, Tukey comparison intervals are recommended as they provide visual consistency with the statistical test being performed and avoid the common confusion where traditional confidence intervals don't overlap but groups share the same significance letter.

Value

An object of class mct (a list with class attributes) containing:

Supported model types

The comparison functions (multiple_comparisons(), pairwise_comparisons() and reference_comparisons()) work with any model for which a get_predictions() method is defined. These are currently:

ARTool (art) models are supported by resplot() but not by the comparison functions: the aligned rank transform makes mean-based comparisons inappropriate. Use ARTool::art.con() for contrasts on ART models instead.

sommer mmer models (the legacy interface) are supported by resplot() but not by the comparison functions: current sommer provides no predict() method for mmer. Refit with sommer::mmes() to use the comparison functions.

To add a new engine, write a ⁠get_predictions.<class>()⁠ method returning a list with elements predictions, sed, df, ylab and aliased_names (plus emmeans_grid for engines backed by emmeans::emmeans()), and add a row to the table above.

References

Jørgensen, E. & Pedersen, A. R. (1997). How to Obtain Those Nasty Standard Errors From Transformed Data - and Why They Should Not Be Used.

See Also

pairwise_comparisons() for testing a chosen subset of pairwise differences as a tidy table, or reference_comparisons() for testing treatments against a chosen reference or control. For guidance on choosing between the two and on multiplicity adjustments, see vignette("choosing-multiple-comparisons", "biometryassist").

Examples

# Fit aov model
model <- aov(Petal.Length ~ Species, data = iris)

# Display the ANOVA table for the model
anova(model)

# Determine ranking and groups according to Tukey's Test (with Tukey intervals)
pred.out <- multiple_comparisons(model, classify = "Species")

# Display the predicted values table
pred.out

# Access the p-value matrix
pred.out$pairwise_pvalues

# Access the HSD value
pred.out$hsd

# Show the predicted values plot
autoplot(pred.out, label_height = 0.5)

# Use traditional confidence intervals instead of Tukey comparison intervals
pred.out.ci <- multiple_comparisons(model, classify = "Species", int.type = "ci")
pred.out.ci

# Plot without confidence intervals
pred.out.none <- multiple_comparisons(model, classify = "Species", int.type = "none")
autoplot(pred.out.none)

# Use a different p-value adjustment instead of Tukey's HSD
multiple_comparisons(model, classify = "Species", adjust = "fdr")

# `by`: run the comparisons independently within each level of another factor.
# Here a 2 x 3 factorial - compare tension levels within each wool type.
m_wb <- aov(breaks ~ wool * tension, data = warpbreaks)
mc_by <- multiple_comparisons(m_wb, classify = "wool:tension", by = "wool")
mc_by
autoplot(mc_by) # faceted by wool

# AOV model example with transformation
my_iris <- iris
my_iris$Petal.Length <- exp(my_iris$Petal.Length) # Create exponential response
exp_model <- aov(Petal.Length ~ Species, data = my_iris)

resplot(exp_model) # Residual plot shows problems

# Fit a new model using a log transformation of the response
log_model <- aov(log(Petal.Length) ~ Species, data = my_iris)

resplot(log_model) # Looks much better

# Display the ANOVA table for the model
anova(log_model)

# Back transform, because the "original" data was exponential
pred.out <- multiple_comparisons(log_model, classify = "Species",
                                   trans = "log")

# Display the predicted values table
pred.out

# Show the predicted values plot
autoplot(pred.out, label_height = 15)

## Not run: 
# ASReml-R Example
library(asreml)

# Fit ASReml-R model
model.asr <- asreml(yield ~ Nitrogen + Variety + Nitrogen:Variety,
                    random = ~ Blocks + Blocks:Wplots,
                    residual = ~ units,
                    data = asreml::oats)

wald(model.asr) # Nitrogen main effect significant

#Determine ranking and groups according to Tukey's Test
pred.out <- multiple_comparisons(model.obj = model.asr, classify = "Nitrogen",
                    descending = TRUE)

print(pred.out, decimals = 5)

# Example using a box-cox transformation
set.seed(42) # See the seed for reproducibility
resp <- rnorm(n = 50, 5, 1)^3
trt <- as.factor(sample(rep(LETTERS[1:10], 5), 50))
block <- as.factor(rep(1:5, each = 10))
ex_data <- data.frame(resp, trt, block)

# Change one treatment random values to get significant difference
ex_data$resp[ex_data$trt=="A"] <- rnorm(n = 5, 7, 1)^3

model.asr <- asreml(resp ~ trt,
                    random = ~ block,
                    residual = ~ units,
                    data = ex_data)

resplot(model.asr)

# lambda = 1/3 from MASS::boxcox()
model.asr <- asreml(resp^(1/3) ~ trt,
                    random = ~ block,
                    residual = ~ units,
                    data = ex_data)

resplot(model.asr) # Look much better

pred.out <- multiple_comparisons(model.obj = model.asr,
                                 classify = "trt",
                                 trans = "power", power = (1/3))

pred.out
autoplot(pred.out, label_height = 0.5)

## End(Not run)

Count Unique Values

Description

Internal helper to count the number of distinct values in a vector. Works for numeric, character, and factor vectors.

Usage

n_unique(x, na.rm = FALSE)

Arguments

Value

Integer. The number of unique values in x (including NA as one distinct value when na.rm = FALSE).

Compare installed version of ASReml-R with available versions

Description

Compare installed version of ASReml-R with available versions

Usage

newer_version(manifest = fetch_manifest(), warn = TRUE)

Value

TRUE if a newer version is available online, FALSE otherwise

Note when per-comparison CIs and adjusted p-values can disagree

Description

Shared by pairwise_comparisons() and reference_comparisons(). For non-simultaneous adjustments the confidence intervals are per-comparison (at level sig) while the p-values are adjusted for multiplicity, so an interval can exclude zero while the adjusted p-value is not significant (or, more rarely, the reverse). Emits a one-time explanatory message() when this actually occurs in the result. Dunnett intervals are the simultaneous intervals and agree with the test by construction, so they are never flagged; "none" likewise never disagrees.

Usage

note_ci_padjust_mismatch(x, sig, method)

Arguments

Value

Invisibly TRUE if a message was shown, otherwise FALSE.

Build predictions, SED and df from an emmeans-backed model

Description

Shared core for the emmeans-backed get_predictions() methods (aovlist, afex_aov, ...). Given the emmeans reference grid for classify, it builds the predicted means, the comparison-specific (matrix) SED and degrees of freedom from the pairwise contrasts, and processes aliased levels. The terms check and ylab are computed by the caller (these differ per engine) and passed in.

Usage

predictions_from_emmeans(model.obj, classify, model_terms, ylab)

Arguments

Value

A list with elements predictions, sed, df, ylab, aliased_names and emmeans_grid.

Print output of multiple_comparisons

Description

Print output of multiple_comparisons

Usage

## S3 method for class 'mct'
print(x, decimals = 2, ...)

Arguments

Value

The original object invisibly.

See Also

multiple_comparisons()

Examples

dat.aov <- aov(Petal.Width ~ Species, data = iris)
output <- multiple_comparisons(dat.aov, classify = "Species")
print(output)
print(output, decimals = 4)

Process aliased treatments in predictions

Description

Process aliased treatments in predictions

Usage

process_aliased(
  pp,
  sed,
  classify,
  exclude_cols = c("predicted.value", "std.error", "df", "Names"),
  vcov = NULL
)

Arguments

Value

List containing processed predictions, sed matrix, aliased names and (when supplied) the subset vcov.

Function to suppress output if desired, especially useful for ASReml output

Description

Function to suppress output if desired, especially useful for ASReml output

Usage

quiet(x)

Arguments

Value

The invisible output of the function called.

Remove existing ASReml installation

Description

Remove existing ASReml installation

Usage

remove_existing_asreml(verbose = FALSE)

Arguments

Produce residual plots of linear models

Description

Produces plots of residuals for assumption checking of linear (mixed) models.

Usage

resplot(
  model.obj,
  shapiro = TRUE,
  call = FALSE,
  label.size = 10,
  axes.size = 10,
  call.size = 9,
  onepage = FALSE,
  onepage_cols = 3,
  mod.obj
)

Arguments

Value

A ggplot2 object containing the diagnostic plots.

Supported model types

resplot() produces residual diagnostics for any model with an extract_model_info() method. These are currently:

This set differs slightly from the comparison functions (see get_predictions()): resplot() additionally supports ARTool (art) models and sommer's legacy mmer interface. Neither is available for the comparison functions — ART uses aligned ranks (use ARTool::art.con()), and current sommer provides no predict() for mmer (refit with sommer::mmes()).

To add a new engine, write an ⁠extract_model_info.<class>()⁠ method returning a list with elements facet, facet_name, resids, fits, k and model_call, and add a row to the table above.

Examples

dat.aov <- aov(Petal.Length ~ Petal.Width, data = iris)
resplot(dat.aov)
resplot(dat.aov, call = TRUE)

Residual plots of linear models.

Description

Produces plots of residuals for assumption checking of linear (mixed) models.

Usage

resplt(model.obj, shapiro = TRUE, call = FALSE, label.size = 10,
axes.size = 10, call.size = 9, mod.obj)

Arguments

Value

A list containing ggplot2 objects which are diagnostic plots.

See Also

biometryassist-deprecated

Produces a skeletal ANOVA table

Description

Produces a skeletal ANOVA table

Usage

satab(design.obj)

Arguments

Value

Prints skeletal ANOVA table to console output.

Setup Colour Palette for Plotting

Description

Internal helper function to generate or validate colour palettes for experimental design plots. Supports predefined palettes (ColorBrewer, Viridis) or custom colours.

Usage

setup_colour_palette(palette, n)

Arguments

Value

Character vector of hex colour codes of length n

Shapiro-Wilk test

Description

Shapiro-Wilk test

Usage

shapiro_test(group_residuals)

Arguments

Visualise a graphical summary of variables from a data frame

Description

Variables are plotted in different ways according to the number of explanatory variables provided as input.

Usage

summary_graph(data, response, exp_var, resp_units = "")

Arguments

Details

With a single explanatory variable, a boxplot grouped by exp_var is produced. With two explanatory variables, a dot-plot with lines connecting the mean of each group is produced, with the first element of exp_var used as the x axis variable, and the second is used to colour the points. Three explanatory variables produces the same as two, but with the third used to facet the plot.

Value

A ggplot2 plot object

Examples

summary_graph(iris, "Petal.Length", "Species", "mm")

# Multiple explanatory variables can be provided as a vector
summary_graph(npk, "yield", c("N", "P"), "lb/plot")

summary_graph(npk, "yield", c("N", "P", "K"), "lb/plot")

Use biometryassist analysis templates

Description

use_template() copies a pre-built analysis template from the biometryassist package to your working directory and optionally opens it for editing. These templates provide standardized approaches for common agronomic analyses.

Usage

use_template(
  template_name = "mixed_model_template.R",
  dest_dir = ".",
  open = TRUE,
  overwrite = FALSE,
  output_name = NULL
)

Arguments

Details

This function is designed to help users get started with biometryassist analyses by providing tested, documented templates. The templates include:

• Suggested package loading

• Commented code explaining steps

• Example data exploration

• Common analysis workflows

If a file with the same name already exists in the destination directory, the function will not overwrite it unless overwrite = TRUE is specified. In this case, it will still open the existing file if open = TRUE.

Value

The file path to the copied template (invisibly). Called primarily for its side effects of copying and optionally opening the template file.

See Also

• list_templates() to see available templates

Examples

## Not run: 
# Copy and open the default analysis template
use_template()

# Copy a specific template without opening
use_template("anova_template.R", open = FALSE)

# Copy to a specific directory
use_template("mixed_model_template.R", dest_dir = "analyses")

# Overwrite an existing file
use_template("mixed_model_template.R", overwrite = TRUE)

## End(Not run)

Calculate the variogram data frame for a model

Description

Calculate the variogram data frame for a model

Usage

vario_df(model.obj, Row = NA, Column = NA)

Arguments

Value

A data frame with the variogram for a model. The data frame contains the spatial coordinates (typically row and column), the gamma for that position and the number of points with the separation.

Examples

## Not run: 
library(asreml)
oats <- asreml::oats
oats <- oats[order(oats$Row, oats$Column),]
model.asr <- asreml(yield ~ Nitrogen + Variety + Nitrogen:Variety,
                    random = ~ Blocks + Blocks:Wplots,
                    residual = ~ ar1(Row):ar1(Column),
                    data = oats)
vario_df(model.asr)

## End(Not run)

Build the 2D variogram heatmap

Description

Internal helper for variogram(). Builds the deterministic 2D heatmap/contour ggplot2 panel from an interpolated grid (see vario_interp()). This is the lower panel of the composite variogram plot.

Usage

vario_ggplot(gdat, row, column, palette)

Arguments

Value

A ggplot2 object.

Interpolate variogram points onto a regular grid

Description

Internal helper for variogram(). Takes the variogram points for a single group (as produced by vario_df()) and interpolates the gamma surface onto a regular n by n grid using pracma::interp2(). The returned grid is the shared input to both the 2D heatmap (vario_ggplot()) and the 3D wireframe.

Usage

vario_interp(points, n = 40)

Arguments

Value

A data frame with columns x, y and z (the interpolated gamma), containing n * n rows.

Display variogram plots for spatial models

Description

Produces variogram plots for checking spatial trends.

Usage

variogram(
  model.obj,
  row = NA,
  column = NA,
  horizontal = TRUE,
  palette = "rainbow",
  onepage = FALSE
)

Arguments

Value

A ggplot2 object.

References

S. P. Kaluzny, S. C. Vega, T. P. Cardoso, A. A. Shelly, "S+SpatialStats: User’s Manual for Windows® and UNIX®" Springer New York, 2013, p. 68, https://books.google.com.au/books?id=iADkBwvario_pointsQBAJ.

A. R. Gilmour, B. R. Cullis, A. P. Verbyla, "Accounting for Natural and Extraneous Variation in the Analysis of Field Experiments." Journal of Agricultural, Biological, and Environmental Statistics 2, no. 3, 1997, pp. 269–93, https://doi.org/10.2307/1400446.

Examples

## Not run: 
library(asreml)
oats <- asreml::oats
oats <- oats[order(oats$Row, oats$Column),]
model.asr <- asreml(yield ~ Nitrogen + Variety + Nitrogen:Variety,
                    random = ~ Blocks + Blocks:Wplots,
                    residual = ~ ar1(Row):ar1(Column),
                    data = oats)
variogram(model.asr)

## End(Not run)