Skip to contents

Computes a pairwise comparison criterion (HSD, LSD, or Bonferroni-corrected LSD) for predicted values obtained from an ASReml-R model, optionally within groups defined by a subset of the factors in term.

For each group the function computes the average SED across all \(\binom{n_g}{2}\) pairs of predictions:

$$ \overline{\text{SED}} = \frac{1}{\binom{n_g}{2}} \sum_{i < j} \sqrt{V_{ii} + V_{jj} - 2V_{ij}} $$

where \(V_{ij}\) is the \((i,j)\) element of the prediction error variance-covariance matrix. The comparison criterion is then:

"HSD" (Tukey)

\(\displaystyle\frac{\overline{\text{SED}}}{\sqrt{2}} \times q_{1-\alpha}(n_g,\, \nu)\), where \(q\) is the Studentised range quantile with \(n_g\) means and \(\nu = \texttt{model\$nedf}\) error df. Two means are declared significantly different if their absolute difference exceeds the HSD.

"LSD" (unadjusted)

\(\overline{\text{SED}} \times t_{\alpha/2,\,\nu}\). No multiplicity correction; liberal for large numbers of comparisons.

"Bonferroni"

\(\overline{\text{SED}} \times t_{\alpha/(2m),\,\nu}\), where \(m = \binom{n_g}{2}\) is the number of pairs. Controls the family-wise error rate.

Usage

compare(
  model,
  term,
  by = NULL,
  type = c("HSD", "LSD", "Bonferroni"),
  pev = TRUE,
  alpha = 0.05,
  ...
)

Arguments

model

An ASReml-R V4 model object.

term

Character string specifying the classify term, e.g. "Treatment:Genotype". Predictions are obtained via asreml::predict.asreml().

by

Grouping specification. Comparisons are performed independently within each level of the group. Accepted forms:

NULL (default)

All predictions form one group.

Single string "Site"

Split by the Site factor in term.

Colon string "Site:Treatment"

Combine Site and Treatment into a single grouping variable and split by their interaction.

Character vector c("Site", "Treatment")

Equivalent to the colon form above.

All by variables must be present in term. At least one variable in term must remain for comparisons to be made.

type

Comparison type. One of "HSD" (Tukey, default), "LSD" (unadjusted), or "Bonferroni".

pev

Logical. If TRUE (default) the prediction error variance (PEV) from predict()$vcov is used. If FALSE and term is in the random part of the model, the posterior variance \(\bm{G} - \text{PEV}\) is used instead.

alpha

Significance level. Default 0.05.

...

Additional arguments forwarded to asreml::predict.asreml().

Value

A data frame with one row per non-missing predicted value, containing:

Factor columns

All factor variables from term.

predicted.value

The BLUE or BLUP.

std.error

Prediction standard error.

{type}

The comparison criterion (constant within each group). Two predicted values in the same group are significantly different if their absolute difference exceeds this value.

avsed

The average SED for the group (constant within group).

See also

waldTest() for pairwise p-values, asreml::predict.asreml()

Examples

if (FALSE) { # \dontrun{
## HSD for Treatment:Genotype, one group per Site
res <- compare(model,
               term = "Treatment:Site:Genotype",
               by   = "Site",
               type = "HSD")

## LSD computed within each Treatment x Site combination
res2 <- compare(model,
                term = "Treatment:Site:Genotype",
                by   = c("Treatment", "Site"),
                type = "LSD")

## Equivalently using colon notation
res3 <- compare(model,
                term = "Treatment:Site:Genotype",
                by   = "Treatment:Site",
                type = "LSD")

## Bonferroni-corrected LSD, single global group
res4 <- compare(model,
                term  = "Treatment:Genotype",
                type  = "Bonferroni",
                alpha = 0.05)

## Posterior variance (BLUPs)
res5 <- compare(model,
                term = "Treatment:Genotype",
                by   = "Treatment",
                pev  = FALSE)
} # }