Skip to contents

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


  sig = 0.05,
  int.type = "ci",
  trans = NA,
  offset = NA,
  power = NA,
  decimals = 2,
  descending = FALSE,
  plot = FALSE,
  label_height = 0.1,
  rotation = 0,
  save = FALSE,
  savename = "predicted_values",



An ASReml-R or aov model object. Will likely also work with lme (nlme::lme()), lmerMod (lme4::lmer()) models as well.


Name of predictor variable as string.


The significance level, numeric between 0 and 1. Default is 0.05.


The type of confidence interval to calculate. One of ci, 1se or 2se. Default is ci.


Transformation that was applied to the response variable. One of log, sqrt, logit, power or inverse. Default is NA.


Numeric offset applied to response variable prior to transformation. Default is NA. Use 0 if no offset was applied to the transformed data. See Details for more information.


Numeric power applied to response variable with power transformation. Default is NA. See Details for more information.


Controls rounding of decimal places in output. Default is 2 decimal places.


Logical (default FALSE). Order of the output sorted by the predicted value. If TRUE, largest will be first, through to smallest last.


Automatically produce a plot of the output of the multiple comparison test? Default is FALSE. This is maintained for backwards compatibility, but the preferred method now is to use autoplot(<multiple_comparisons output>). See autoplot.mct() for more details.


Height of the text labels above the upper error bar on the plot. Default is 0.1 (10%) of the difference between upper and lower error bars above the top error bar.


Rotate the text output as Treatments within the plot. Allows for easier reading of long treatment labels. Number between 0 and 360 (inclusive) - default 0


Logical (default FALSE). Save the predicted values to a csv file?


A file name for the predicted values to be saved to. Default is predicted_values.


Deprecated. Use descending instead.


Deprecated. Predicted values are calculated within the function from version 1.0.1 onwards.


Deprecated. Use classify instead.


Other arguments passed through to predict.asreml().


A list containing a data frame with predicted means, standard errors, confidence interval upper and lower bounds, and significant group allocations (named predicted_values), as well as a plot visually displaying the predicted values (named predicted_plot). If some of the predicted values are aliased, a warning is printed, and the aliased treatment levels are returned in the output (named aliased).


Some transformations require that data has a small offset applied, otherwise it will cause errors (for example taking a log of 0, or square root of negative values). In order to correctly reverse this offset, if the trans argument is supplied, an offset value must also be supplied. If there was no offset required for a transformation, then use a value of 0 for the offset argument.


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


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

# Display the ANOVA table for the model
#> Analysis of Variance Table
#> Response: Petal.Length
#>            Df Sum Sq Mean Sq F value    Pr(>F)    
#> Species     2 437.10 218.551  1180.2 < 2.2e-16 ***
#> Residuals 147  27.22   0.185                      
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

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

# Display the predicted values table
#>      Species predicted.value std.error  df groups  ci  low   up
#> 1     setosa            1.46      0.06 147      a 0.1 1.36 1.56
#> 2 versicolor            4.26      0.06 147      b 0.1 4.16 4.36
#> 3  virginica            5.55      0.06 147      c 0.1 5.45 5.65

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

if (FALSE) {
# ASReml-R Example

#Fit ASReml 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, 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)


# Perform Box-Cox transformation and get maximum value
out <- MASS::boxcox(ex_data$resp~ex_data$trt)
out$x[which.max(out$y)] # 0.3838

# Fit cube root to the data
model.asr <- asreml(resp^(1/3) ~ trt,
                    random = ~ block,
                    residual = ~ units,
                    data = ex_data)
resplot(model.asr) # residual plots look much better

#Determine ranking and groups according to Tukey's Test
pred.out <- multiple_comparisons(model.obj = model.asr,
                                 classify = "trt",
                                 trans = "power", power = (1/3))