multiple_comparisons {biometryassist}

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 = NA,
  offset = NA,
  power = NA,
  decimals = 2,
  descending = FALSE,
  plot = FALSE,
  label_height = 0.1,
  rotation = 0,
  save = FALSE,
  savename = "predicted_values",
  order,
  pred.obj,
  pred,
  ...
)

Arguments

model.obj

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

classify

Name of predictor variable as string.

sig

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

int.type

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

trans

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

offset

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.

power

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

decimals

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

descending

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

plot

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.

label_height

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.

rotation

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

save

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

savename

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

order

Deprecated. Use descending instead.

pred.obj

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

pred

Deprecated. Use classify instead.

...

Other arguments passed through to predict.asreml().

Details

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.

Value

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

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. https://pure.au.dk/portal/en/publications/how-to-obtain-those-nasty-standard-errors-from-transformed-data–and-why-they-should-not-be-used(d649ca20-d15f-11db-8e26-000ea68e967b).html

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
pred.out <- multiple_comparisons(model, classify = "Species")

# Display the predicted values table
pred.out

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



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

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

pred.out

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

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

pred.out
autoplot(pred.out)

## End(Not run)